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ABSTRACT 

Cebeci's  viscous  inviscid  interaction  program  was  applied  to  the  analysis  of 
steady  two  dimensional  incompressible  flow  past  a  NACA  65-213  airfoil  at  zero  angle 
of  attack  at  a  Reynolds  number  of  240,000.  Predicted  boundary  layer  characteristics 
were  found  to  be  quite  sensitive  to  the  choice  of  boundary  layer  transition  begin  and 
length.  Good  agreement  with  the  experimental  results  of  Hoheisel  er  al  could  be 
obtained  by  proper  choice  of  both  transition  begin  and  length. 
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I.  INTRODUCTION 

The  prediction  of  the  stall  characteristics  is  of  major  importance  in  aeronautical 
engineering  to  determine  the  operating  limits  of  an  aircraft.  Therefore,  the 
development  of  reliable  and  accurate  numerical  methods  for  predicting  separated  flow 
regions  is  one  of  the  most  challenging  problems  in  Computational  Fluid  Dynamics 
(CFD). 

In  this  thesis  we  limit  ourselves  to  the  problem  of  incompressible  two- 
dimensional  airfoil  flows.  Two  basic  methods  are  available  to  compute  viscous  flows 
which  include  regions  of  How  separation.  The  first  approach  is  based  upon  a  solution 
of  the  full  Navier-Stokes  equations  (or  some  approximate  form,  such  as  the  parabolized 
Navier-Stokes  equations).  This  approach  has  the  disadvantage  of  being  very  expensive 
and  time  consuming.  The  second  approach  is  based  upon  the  so-called  Viscous-Inviscid 
Interaction  Method.  The  outer  flow  is  computed  using  the  inviscid  flow  equations.  The 
inner  flow  (close  to  the  airfoil)  is  obtained  from  a  numerical  solution  of  the  Prandtls 
boundary  layer  equation.  However,  in  contrast  to  the  well-known  classical  boundary 
layer  computations  the  pressure  cannot  be  prescribed  a  priori,  but  must  be  found 
iteratively  (i.e.,  by  viscous-inviscid  interaction).  This  approach  has  the  advantage  of 
being  much  faster  and  more  efficient  than  the  Navier-Stokes  solutions. 

The  approach  chosen  in  this  thesis  is  based  upon  the  viscous-inviscid  interaction 
method  developed  by  T.  Cebeci  and  collaborators  at  the  Douglas  Aircraft  Company.  In 
chapter  2  the  fundamental  equations  are  summarized  .  Chapter  3  is  devoted  to  a 
discussion  of  the  inviscid  flow  method,  especially  the  so-called  Panel  Method  first 
introduced  [Ref.  4]  by  Hess  and  Smith.  In  chapter  4  the  direct  and  interactive 
boundary  layer  methods  are  discussed,  followed  by  a  brief  explanation  of  the 
turbulence  model  used.  Finally,  the  computer  program  is  explained  and  computed 
results  are  presented  for  the  NACA  65-213  airfoil  and  compared  with  detailed 
measurements  by  Hoheisel  et  al.    [Ref.  14]. 
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-     II.  FUNDAMENTAL  EQUATIONS 

This  chapter  presents  the  equations  used  in  our  analysis,  which  involve  the 
development  of: 

I.     Conservation  of  mass  (Continuity  equation* 

Conservation  of  momentum  (Newton's  Second  Law) 

Inviscid  How  equations 

Boundary  layer  equation 

Turbulence  model 

With  these  equations  we  can  predict  the  behavior  of  a  body  moving  through  the  fluid. 
In  our  case  we  -ire  dealing  with  an  airfoil  in  two  dimensional,  steady,  inviscid  and 
viscous  flow. 


A.       CONSERVATION  OF  MASS 

Let  us  apply  the  principle  of  conservation  of  mass  to  a  small  volume  of  space 
through  which  the  fluid  can  move  freely.  For  convenience,  we  shall  use  a  cartesian 
coordinate  system  Cx\y.r).  Furthermore,  in  the  interest  of  simplicity,  we  shall  treat  a 
•2-D  flow,  that  is.  one  in  which  there  is  no  How  along  the  r  -axis.  Flow  patterns  are  the 
same  for  any  x-y  plane.  As  indicated  in  the  sketch  of  Figure  2.1.  the  component  of  the 
fluid  velocity  in  the  x  direction  will  be  designated  by  u,  and  in  the  y  direction  by  v.  The. 
net  outflow  of  mass  through  the  surface  surrounding  the  volume  must  be  equal  to  the 
decrease  of  mass  within  the  volume.  The  mass  ilow  rate  through  the  surface  bounding 
the  element  is  equal  to  the  product  of  the  density,  the  velocity  component  normal  to 
the  surface,  and  area  of  that  surface. 

A  first-order  Taylor  scries  expansion  is  used  to  evaluate  the  flow  properties  at  the  faces 
of  the  element  [Ref.  3]  since  the  properties  are  a  function  of  position.  Consider  the 
How  out  of  the  volume  as  positive,  then  the  net  outflow  of  mass  per-unit  time  is  the 
summation  of 
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Figure  2.1     Sketch  illustrating  the  velocitv 

and  the  densitv  for  mass  flow  balance 

throueh  a  fixed  volume  in  2-D. 
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which  must  equal  the  rare  at  which  the  mass  contained  within  the  element  decreases 

-3i^y  <Z2) 

c/t 


Combining  euuation  2.1  and  2.2.   dividing  bv  A.rAv.  one  sets 


dp        d  Q 


In  vector  form,  ecua'tion  2.3  is 


-— +r.(pl  )  =0  (2.4) 

at 


Because  the  pressure  variations  that  occur  in  relatively  low  speed  flow  are  sufficiently 
small,  the  density  is  essentially  constant.  For  these  incompressibe  flows,  the  continuity 
equation  becomes 


dn       d'' 
ax      dy 

or.  in  vector  form  V  •  T  =  0 


B.       CONSERVATION  OF  MOMENTUM 

The  equation  of  the  conservation  of  linear  momentum  is  obtained  by  applying 
Newton's  2nd  Law  [Ref.  3]  where  the  net  force  acting  on  a  fluid  particle  is  equal  to  the 
time  rate  of  change  of  the  linear  momentum  of  the  fluid  particle.  As  the  fluid  moves  in 
space,  its  shape  and  volume  may  change,  but  its  mass  is  conserved.  Thus,  using  a 
coordinate  system  that  is  neither  accelerating  nor  rotating,  called  an  inertial  coordinate 
svstem.  we  mav  write 
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r         Dv 

Dt 


The  velocity  V  of  cTfluid  particle  is.  in  general,  an  explicit  function  of  time  t  as  well  as 
of  its  position  xy.  Furthermore,  the  position  coordinates  jcvv  of  the  fluid  particle  are 
themselves  a  function  of  time.  Since  the  time  differentiation  of  equation  2."  follows  a 
given  particle  in  its  motion,  the  derivative  is  frequently  termed  the  particle  or 
substantial  derivative  of  V,  since  \\xy.i)  and  x{n,y{i). 


—  u f-  v 1 ( 2.S ) 


Dt  dx  dy        dt 

where,  u  =  ~   -     and         v  =  - 
dt  at 

Therefore,  the  acceleration  of  a  fluid  particle  is 


DV       dV         . 


Thus,  the  substantial  derivative  is  the  sum  of  the  local,  time  dependent  changes 
that  occur  at  a  point  in  the  flow  field  and  of  the  convection  in  space.  When  the  local. 
time  dependent  changes  are  zero.  dVcr=0  ,  such  flows  are  known  as  steady-state 
flows.  The  principal  forces  that  act  on  the  body  are  those  which  act  directly  on  the 
mass  of  the  fluid  element,  the  body  forces  ,  and  those  which  act  on  its  surface,  the 
pressure  forces  ami  shear  forces  known  as  surface  forces.  The  stress  system  acting  on  an 
element  o[  the  surface  is  illustrated  in  Figure  2-2.  The  stress  components  t  acting  on 
the  small  cube  are  assigned  subscripts.  The  first  subscript  indicates  the  direction  of  the 
normal  direction  to  the  surface  on  which  the  stress  acts  and  the  second  subscript 
indicates  the  direction  in  which  the  stress  acts.  Thus.  T  denotes  a  stress  acting  in  the 
y  direction  on  the  surface  whose  normal  points  in  the  \  direction.  The  properties  of 
most  fluids  have  no  preferred  direction  in  space:  that  is.  fluids  are  inotropic. 
Again,  the  stresses  on  the  fluid  element  can  be  obtained  from  a  Taylor  scries 
expansion.     In  general,  the  various  stresses  change  from  point  to  point.    Thus,  they 
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Fisure  2.2    Stresses  acting  on  a  2-D  element  of  fluid. 
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produce  net  forces  on  the  fluid  particle,  which  cause  it  to  accelerate.  The  forces  acting 
on  each  surface  are  obtained  by  taking  into  account  the  variations  of  stress  with 
position,  by  using  t!Te  center  of  the  element  as  a  reference  point. 

To  simplify  the  illustration  of  the  force  balance  on  the  fluid  particle  we  shall 
again  consider  a  2-D  flow,  as  indicated  in  Figure  2-2.  The  resultant  force  in  the  x- 
dircction.  for  one  unit  lencth  in  z  is 


dx  dy 


(2.10) 


Where  /'  is  the  bodv  force  per-unit  mass  in  the  x  direction.  The  most  common  bodv 
force  for  the  flow  fields  is  that  o[  gravity.  Equation  2.10  is  the  left  hand  side  of 
equation  2.".  For  the  right  hand  side,  combine  mass  term  with  equation  2.10  in  the  x 
direction 


Du      ,  ■ 
m—  =  [pAxAy) 


at 


(•2.11) 


From  equation  2.10  and  2:1 1.  substitute  into  equation  2.7  divided  by  Aa*Av  we  have 


-         9  d  \  dn 


dx  dy 


dt 


similarlv  in  the  v  direction 


dx 


d_ 


fd> 


-'  =  '\--{V 


(2.13) 


Next,  we  need  a  relation  between  the  stresses  and  the  motion  of  the  fluid.  For  a 
fluid  at  rest  or  for  inviscid  fluid  motion,  there  is  no  shearing  stress  and  the  normal 
stress  is  in  the  nature  of  a  pressure.  Tor  fluid  particles,  the  stress  is  related  to  the  rate 
of  strain  by  a  physical  law  based  on  the  following  assumptions: 
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1.  Stress  components  may  be  expressed  as  a  linear  function  of  the  components  of 
the  rate  of  the  strain.  The  friction  law  for  1-D  flow  of  a  Newtonian  fluid  is  a 
special  case  of  this  linear  stress/rate  of  strain  relation,  i.e.,  T  =  ]i  (  cu.dy  ), 
where  \i  is  the  fluid  viscosity. 

2.  The  relation  between  the  stress  components  and  strain  rate  components  must 
be  invariant  to  a  coordinate  transformation  consisting  of  either  a  rotation  or  a 
mirror  reflection  of  axes,  since  a  physical  law  cannot  depend  upon  the  choice  of 
the  coordinate  system. 

3.  When  all  velocity  gradients  are  zero  (i.e.,  the  shear  stress,  vanishes),  the  stress 
components  must  reduce  to  the  hydrostatic  pressure  p. 

For  a  fluid  that  satisfies  these  criteria,  in  two  dimensional  flow 


dn       2 
ox       3 


2.14) 


dU 


2.15] 


rxy  =  v 


dn       dr 
.  dy      dx 


2.1C) 


with  the  appropriate  expressions  for  the  surface  stresses,  substitute  into  equation  2.12 
and  2.13.  we  have 


dn 

—  +  (V-V)n 

dt 


'1l 

dt 


+  (V-V! 


=  pA- 

dp       d 

dx  +  dx 

'      dn 

2a<- 

ox 

0 

3M 

•V) 

+ 

d 

,  dn       dr   ' 

du 

dy       dx 

=  Pfy  ~ 

dJL+L 
dy  '   dx 

,  dn       dr . 
.      dy      dx   . 

+ 

d 
dy 

V 

h 

h 

2 
T  3 

■/'(v-r) 

(2.17 


(2.18) 
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These  general  differential   equations  for  the  conservation   of  linear  momentum  are 
known  as  the  Navier  -  Stokes  equations.    When  we  are  dealing  with  incompressible  and 
2-D  flow,  then  fronTequation  2.6.  V.T  =  0,  and  the  body  forces  arc  neglected,  and  r 
+  r,.,.  =  -  2p.   Therefore  equation  2.17  and  2. IS  can  be  written  as 


dx 


dv 
ox 


dn       du 

1  dp 

d'u       d'u 

du       dt 

p  ox 

,  dx2       dij' 

dv      dv 
dij      dt  " 

1  dp 

-  +  u 

P   $U 

r  d2v       d2v 
.  dx2       dij2 

(2.10) 
(2.20) 


These  are  the  well  known  Xavier  •  Stokes  equation  for  2-D,  incompressible,  viscous 


flow. 


C.       INVISCID  FLOW  EQUATION 

Inviscid  flow  represents  an  ideal  flow,  where  the  effects  of  viscosity  are  zero.  In 
reality,  this  is  not  true  because  every  medium  has  viscosity,  eventhough  it  may  be  very 
small. 

Why  is  the  inviscid  flow  important  in  dealing  with  a  body  moving  in  the  viscous 
fluid  ? 

In  1904  L.Prandtl  came  up  with  the  answer.  For  high  Reynolds  Number  on  a 
body  moving  in  a  viscous  fluid,  two  regions  can  be  distinguished.  The  effect  of 
viscosity  can  be  neglected  outside  a  very  thin  region  near  the  body  which  is  called 
boundary  layer.  Inside  the  boundary  layer  the  viscous  effects  are  important  (further 
discussion  in  viscous  flow  section).  This  is  the  reason  why  the  inviscid  flow  remains 
important  in  the  computation  of  fluid  dynamics,  eventhough  it  represents  an  ideal  case. 

1.   Potenial  Flow 

Since  the  flow  upstream  of  the  body  is  uniform  then  it  is  also  irrotational  (  in 
inviscid  flow). 


i  =  T  x  V 

Consider  the  following  vector  identity  ,  if  (p  is  a  scalar  function,  then 
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(2.21) 


fx(V^)  =  0  (122) 

i.e..  the  curl  of  the  gradient  of  a  scalar  function  is  identically  zero.  Comparing  equation 
2.21  and  2.22.  we  see  that 

f  =  V?  (2.23) 

Equation  2.23  states  that  for  an  irrotational  flow  there  exists  a  scalar  function  (p  such 
that  the  velocity  is  given  by  the  gradient  of  (p.  We  denote  (p  as  the  velocity  potential,  (p 
is  a  function  of  the  spatial  coordinates,  i.e..  (p  =  (p(.r,y).  And  from  the  definition  of  the 
gradient  in  cartesian  coordinates  .  equation  2.23  can  be  written  as 


(2.24) 


dj 

dc 

=  — — 

aud 

r  =         * 

dx 

dy 

Thus,  (p  has  the  property  that  its  partial  derivative  in  any  direction  is  the 
velocity  component  in  that  direction.  It  follows  that  the  existence  of  <p  is  the  sole 
criterion  for  irrotationaiity.  The  usefulness  of  the  velocity  potential  in  flows  of 
practical  significance  derives  from  the  circumstance  that,  for  a  body  in  relative  motion 
in  an  originally  irrotational  flow,  the  circulation  vanishes  around  any  contour  that  does 
not  include  the  body  or  does  not  intersect  the  boundary  layer  or  the  wake,  therefore,  a 
velocity  potential  can  be  found  to  describe  the  flow  everywhere  outside  the  boundary 
layer  or  the  wake.  When  a  flow  field  is  irrotational.  hence  allowing  the  velocity 
potential  to  be  defined,  there  is  a  tremendous  simplification.  Instead  of  dealing  with  the 
velocity  components  (say.  u.v  and  w)  as  the  unknowns,  hence  requiring  three  equations 
for  three  unknowns,  we  can  instead  deal  with  the  velocity  potential  as  one  unknown. 
therefore,  requiring  the  solution  of  only  one  equation  for  the  flow  field  and  the  velocity 
components  can  be  obtained  from  equation  2.2-4.  Because  irrotational  flows  can  be 
described  by  the  velocity  potential  (p  such  flows  are  called  potential  flows. 
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2.  Governing  Equation  for  irrotational.  incompressible  flow  :  Laplace  Equation. 

We  have  seen  from  equation  2.6.  that  the  conservation  of  mass  for  an 
incompressible  flcmTakes  the  form  V.V  =  0.  In  addition,  for  irrotational  How  we  have 
seen  in  equation  2.23  ,  V  =  V<p.  Therefore,  for  a  How  that  is  both  incompressible  and 
also  irrotational.  equations  2.6  and  2.23  can  be  combined  to  yield. 


Vv  =  0.         or         —  T  — =  0  U-o) 

di-       dy- 


Equation  2.25  is  called  Laplace's  Equation  .  one  of  the  most  extensively  studied 
equations  in  mathematical  physics. 

Note  that  Laplace's  equation  is  a  second  order  linear  partial  differential 
equation.  The  fact  that  it  is  linear  is  particularly  important,  because  the 'superposition 
of  any  particular  solution  of  a  linear  differential  equation  is  also  a  solution  of  the 

equation.  For  example,  if  q> . .  (p.,.  (p, (p    represent  n  separate  solutions  of  equation 

2.25.  then  the  sum  (p  =  (pj   •     <p^  +  4-  cp^  is  also  a  solution  of  equatiqn  2.25. 

Since  irrotational.  incompressible  flow  is  governed  by  Laplace's  equation  and 
Laplace's  equation  is  linear,  we  conclude  -that  a  complicated  flow  pattern  for  an 
irrotational.  incompressible  flow  can  be  synthesized  by  superposition  of  elementary 
flows  which  are  also  irrotational  and  incompressible.  The  singularity  (or  panel) 
methods  presented  in  the  next  chapter  are  based  on  this  idea. 

D.       BOUNDARY  LAYER  EQUATION 

Up  to  this  point,  we  have  been  dealing  with  flow  outside  the  boundary  layer, 
where  viscous  effects  remain  small.  In  the  region  within  the  boundary  layer,  velocity 
gradients  are  high  even  with  very  small  viscosity.  Therefore,  it  becomes  very  important 
to  deal  with  a  real  fluid.  Before  the  boundary  layer  can  be  analyzed  further,  we  have 
to  know  what  governing  equations  can  be  used  in  the  practical  analysis.  The  objective 
is  to  predict  viscous  flows  by  means  of  the  boundary  layer  method,  instead  of  solving 
the  complete  Navier-Stokes  equations. 

From  the  previous  derivations,  equations  2.5  ,  2.19  and  2.20  are  used  here.  In 
order  to  simplify  these  equations  we  have  to  make  some  assumptions:  two- 
dimensional,  steady,  constant  fluid  properties,  and  no  body  forces.    Another  important 
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assumption  is  that  the  boundary  layer  thickness  is  very  small  compared  to  the  length  of 
the  body  (airfoil  in  this  case).  With  these  assumptions.  L.  Prandtl  in  1904  introduced 
the  "order  of  magnTTuJe"  estimate  into  equation  2.5,  2.19  and  2.20.  In  order  to  do  this, 
all  linear  dimensions  will  be  referenced  to  the  characteristic  length  /  and  all  velocities 
will  be  referenced  to  U,  therefore  /  and  U  can  be  said  to  have  order  of  magnitude  of 
one  t  written  as  0(  1 )  )  and  with  the  above  assumptions,  y  will  correspond  to  the 
boundary  layer  thickness  (6).  Then  the  continuity  equation  2.5  can  be  written  as 

U\         (v 


I 


b 


and 


,u 


Because  6  is  very  small,  from  the  above  relation  we  can  imply  that  v,  the  velocity  in 
the  y  direction  must  be  very  small.  Therefore,  the  continuity  equation  still  holds  in  the 
boundary  layer.    From  the  Navier-Stokes  equation  in  the  .y  direction  (eqn  2.19) 


du  du 


1  d_p_ 
p  dx 


+  v 


d2n 


d'' 


du- 


2.10! 


Introducing  the  order  of  magnitude  in  the  above  equation,  we  have 


u 


7- 
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p_ 
pi 


v 


u_ 

6l 


If  all  above  terms  multiplied  by      (  — -  )       then 


pin 


Ul/u 


r-/<v2 

Ul/u 


If  the  Revnolds  Number  is  hieh.  the  fourth  term 


=  0{:ero) 


Ul/u 

I" /i>~        rery  large 

Ul/u  large 

r-       ( 

therefore     —  =  0[  R, 


i=°Uh 


£  zeroB  0(1) 
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So.  equation  2.19  becomes: 


„  d"  d"  1  dp  d'n 

^       c^        p  di       d>r  (126> 


For  the  y  direction  (  equation  2.20  )  in  terms  of  order  of  magnitude  we  obtain 

U'S\  (U-S\  (     \dp\  (    in\  (    U 

V-—         ,      \   v 


/'-    /  \    1-   J  \    pdijj  \     lz  J  \     bl 


all  above  terms  multiplied  by  (   —  )       then 
(  f>2\  (     1   6    P\  (     U£\  /   [7tS 


r-  j    "    v    p^;2  *  J       \  iztT2J    '   \6iu* 


p)    '    {r-J    '    {    pU*J    '       U"2/^    '■   U</" 


a  ■  a)  •  -*  ■  °(a  ■ » * 


^  L* 


The  only  term  left  is  the  pressure  term,  because  all  other  terms  are  of  higher  order 

dp 

J7j-°  (2.27) 

This  is  very  important  because  it  tells  us  that  the  pressure  across  the  boundary  layer 
remains  constant.  The  triplet  of  equations  (  2.5  ,  2.26  and  2.27  )  and  the  boundary 
condition  of  zero  normal  and  tangential  velocity  are  known  as  the  Prundtl's  boundary 
layer  equations. 

E.       TURBULENT  FLOW  EQUATION 

Since  the  continuity  equation  (2.5)  and  the  Navier-Stokes  equation  (2.19)  make 
no  assumptions  regarding  the  type  of  flow,  they  are  instantaneously  valid  in  both  the 
laminar  and  turbulent  flow  regimes.  However,  it  is  too  difficult  to  deal  with 
instantaneous  properties  in  turbulent  flow.  Therefore  we  introduce  the  time-averaged 


->-> 


properties.    In  our  notation,  prime  denotes  the  fluctuation  quantity,  and  bar  denotes 
the  mean  value. 

u  =  n  +  «' 

v  =  P  +  v' 

p  =  P  +  P 

Introducing  these  properties  into  the  continuity  equation  and  averaging,  we  have 


d  a 

—U-  u>)  -r   r-iP+P')  =  0 
dx  ay 


Carrying  out  the  integration  term  by  term  (details  are  in  [Ref.  17]  ).  yields: 


— (n-riiM  -  -t(r-r')  =  0 


dx' 


d'J 


or 


since         r^fl)  =  — («)=  —(a)         aud         r-{«')   =   tH"')  =  ° 
c>.r  tf.r  dx  dx 


dx 


Therefore  the  continuity  equation  for  turbulent  flow  becomes; 


(2.28) 


A  similar  procedure  can  be  applied  to  the  Navier-Stokes  Equation 


dxx  d\j  p  dx         [  dx1       dy1 


-t-( «'«')-  T-(u'f)  (2.29) 


The  last  two  terms  correspond  to  the  normal  and  shear  stress  terms  respectively,  which 
we  call  the  Reynolds  Stresses  or  Turbulent  Stresses. 


III.  PANEL  METHOD 

A.       INTRODUCTION 

The  thin  airfoil  theory  is  just  what  it  says,  it  applies  only  to  a  thin  airfoil  at  small 
angle  of  attack.  It  is  not  much  used  these  days  for  the  analysis  or  design  of  single 
element  airfoils.  It  does  give  fairly  good  results  for  airfoils  of  12  %  thickness  or  less. 
On  the  other  hand,  the  determination  of  the  aerodynamic  characteristics  of  thick, 
highly  cambered,  slotted  surfaces,  with  single  or  multiple  flaps  and  mutual  interference 
effects  among  wings,  fuselage,  nacelles,  and  so  forth,  requires,  in  general,  the  use  oi' 
numerical  methods. 

In  the  1960's  Mess  and  Smith  at  McDonnell  Douglas  introduced  a  method,  the 
so  called  Panel  Method  [Ref.  4],  as  a  numerical  approach  for  2-D  flows  which  can  be 
extended  to  3-D  potential  flow  problems.  Such  methods  are  called  panel  methods 
because  the  body  surface  is  approximated  by  a  collection  of  panels  .  There  are  a 
number  of  ways  to  set  up  the  panel  method.  To  begin  with,  there  are  choices  even  as 
to  the 'type  of  singularity  used,  sources,  doublets,  vortices  or  a  combination  of  source 
and  vortex  distributions. 

Panel  methods  as  a  numerical  approach  for  predicting  forces  and  moments  acting 
on  the  body  gave  goo'd  agreement  with  the  reliable  published  data.  The  application  of 
the  panel  method  requires  that  the  problem  can  be  formulated  such  that 

1.  the  body  can  be  represented  by  a  closed  polygon  of  a  finite  number  of  elements, 
called  panels  connected  by  nodes. 

2.  the  flow  tangency  condition  is  satisfied  in  the  middle  of  the  panels  (control 
points)  to  avoid  the  inaccuracies  of  thin  airfoil  theory. 

3.  the  singularity  distribution  of  each  element  is  approximated  by  some  kind  of 
analytical  function.  Also,  the  singularities  should  be  distributed  on  the  body 
surface  rather  than  on  the  chord  line  or  any  other  line  within  the  body. 

Before  we  discuss  further  details,  it  is  necessary  to  know  about  the  basic  formulation 
for  source  and  vortex  distribution  as  a  singularity  parameter. 
1.  Single  Source  and  Source  Panel 

A  flowfield  where  all  streamlines  are  straight  lines  emanating  from  a  source 
point,  is  called  a  source  flow.    On  the  other  hand,  when  the  flow  direction  is  inward,  it 


24 


is  called  sink.  flow.    7.T  =   (7  at  every  point  except  at  the  origin  where  7.1"  becomes 
infinite.  The  source  flow  is  irrotational  at  even.-  point.    Its  velocity  potential  is  defined 

as  — 


single  source         ~(x.y)=  ^-     In  r  (}l) 


where  A  is  defined  as  the  source  strength  and  r  is  the  distance  from  the  considered 
point  (.rvr)  to  the  source.  When  we  are  dealing  with  non-lifting  flows  over  a  body,  we 
can  superimpose  elementary  source  Hows  in  order  to  obtain  a  complete  solution.  This 
method  is  called  source  panel  method. 


I 
source  panel         f[x. y)  =  —  /       lurds  ,3.2) 


where-/  is  the  panel  length. 

2.  Single  Vortex  and  Vortex  Panel 

The  How  where  all  streamlines  are  concentric  circles  about  a  given  point  is 
defined  as  single  vortex  flow.  The  velocity  along  any  given  circular  streamline  is 
constant,  but  varies  from  streamline  to  streamline.  Its  velocity  potential  can  be  written 
as 


suisrle  vortex        ^\x,y)  =  - 


r 


o- 


(3.3) 


"•her?     0  =  arctaa]  — 

L-f-  *r 


,    with     [jcv.yv)     as  the  center  of  vortex 

Let  us  imagine  a  straight  line  perpendicular  to  the  page  (Figure  3.1).  This  line  is  a 
straight  vortex  filament  of  strength  T  and  the  flow  induced  in  any  planes  perpendicular 
to  the  vortex  of  strength  T .  i.e..  the  flows  in  the  plane  1  to  the  vortex  filament  at  o 
and  o'  arc  identical  to  each  other  and  are  identical  to  the  How  induced  by  a  point 
vortex  of  strength  T . 
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Figure  3.1     Vortex  Sheet  Representation. 
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Refer  to  Figure  3.1,  imagine  an  infinite  number  of  straight  vortex  filaments 
side  by  side,  where  the  strength  of  each  filament  is  infinitesimally  small.    Define  y  = 
y(s)  as  the  strength,  of  the  vortex  sheet  per  unit  length  along  s,  then  the  strength  of  the 
infinitesimal  portion  ds  of  the  sheet  is  yds  and  the  small  section  of  the  vortex  sheet  of 
strength  yds  induces  an  infinitesimally  small  velocity  dv  at  point  P(x,z) 


,  ds 

di-  =  -  -,  -—     r,  j_  t0     r 

2-  (3.4) 

It  is  sometimes  more  convenient  to  use  the  velocity  potential  (p,  and  the  increment  in 
velocity  potential  d(p  induced  at  P{x,z)  by  elemental  vortex  yds 


</-  =  -—     9 


For  the  entire  vortex  sheet,  from  a  to  b 


j 


(3.5) 


(3.6) 


Equation  3.5  is  useful  in  classical  thin  airfoil  theory  and  equation  3.6  is  important  for 
the  numerical  vortex  oanel  method. 


B.       SOURCE  AND  VORTEX  DISTRIBUTION 

Having   introduced   the   basic   idea    of  the   panel   method,   we   can   use   these 
singularities  either  alone  or  in  combination.    This  method  is  due  to  Hess  and  Smith. 

Thus,  the  potential  may  be  decomposed  in  a  manner  such  that 

(3.7) 
with  (pQ  being  the  potential  of  the  uniform  onset  flow,  and  <ps  and  <pv  the  potentials  due 
to  source  and  vortex  distributions  respectively,  hence 

r        J         1-  (3.S) 
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f    -•(») 


frb  (3.9) 


in  which  the  integrations  are  to  he  performed  over  the  body  surface.  Because  of  the 
superposition  principle,  this  (j>  automatically  satisfies  Laplace's  Equation  (  see  equation 
2.25)  and  the  boundary  condition  at  infinity.  It  will  be  the  solution  we  seek,  if  <T)  s )  and 
7<s)  are  determined  so  as  to  meet  the  boundary  condition  of  flow  tangency  and  the 
Kutta  condition  (to  be  discussed  in  the  next  section). 

Hess  and  Smith  assumed  the  vortex  strength  to  be  constant  over  the  body 
surface  and  the  source  strength  must  vary  over  the  surface.  Since  the  Kutta  condition 
involves  only  the  trailing  edge,  the  vortex  strength  can  be  represented  by  a  single 
number.  Thus,  if  one  distributes  on  or  within  the  body  surface,  vortices  whose  net 
strength  is  the  correct  circulation,  the  problem  is  solved  if  sources  can  be  distributed 
over  the  body  surface  so  as  to  make  the  total  velocity  field  (comprised  oi~  the  onset 
ilow  and  the  velocity  fields  due  to  sources  and  vortices)  tangent  to  the  body  surface, 
regardless  of  how  the  vortices  are  distributed.  However,  the  integrals  in  equations  3.S 
and  3.9  are  hard  to  evaluate,  even  for  simple  forms  -of  the  source  and  vortex  strength. 
unless  the  surface  on  uhich  the  sources  and  vortices  are  distributed,  is-  a  straight  line. 
Thus,  we  select  a  certain  number  of  points  on  the  body  contour,  called  nodes,  and 
connect  the  nodes  with  straight  lines,  which  become  the  panels  of  the  method  (see 
Tigure  3.2  ) 

We  then  distribute  the  sources  and  vortices  on  the  straight  line  panels,  so  that  the 
potential  given  by  equation  3.7  can  be  written  as: 


p  =  To(jrcoscr  +  ysiaa)+    £    /  (  —  hi r  -  ~-<9)  <fc  (3  10) 

In  most  cases,  equation  3.10  still  allows  an  exact  solution  of  the  flow  problem.  The 
exceptional  cases  are  those  in  which  the  sources  and  vortices  must  be  distributed 
exactly  on  the  body  surface  ;  to  be  mathematically  precise,  in  which  the  potential 
cannot  be  continued  anlytically  across  the  body  surface.  By  increasing  the  panel 
density,  the  body  shape  can  be  better  approximated.  This  is  the  only  major 
approximation  of  the  panel  method,  one  that  becomes  more  accurate  as  the  number  of 
panels  increases. 
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Figure  3  2     Discretization:  a).  Actual  airfoil 
o).  Airfoil  alter  discretization. 
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Fisure  3.3     Representation  ofi-th  and  j-th  panels 
in  the  Panel  Method. 
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To  implement  this  method,  we  need  a  nomenclature.  Let  the  i-th  panel  be 
defined  as  the  one  between  the  i-th  and  (/+  l)th  nodes,  and  its  inclination  to  the  x  axis 
be  0j.   Therefore  \b&_  relation  between  midpoints  and  boundary  points 

IX  +  A'm.1) 


X,   = 


Ui  = 


r  +  rt+1 


and  the  computation  of  the  angle  9 


(9,  =  arctau 


lUi  -  r. 

L  AVi  -  A'. 


From  the  geometry  in  Figure  3.3  a  relation  for  angle  and  distance  between  two  panels 
can  be  obtained 


rtJ  =  J[xi  -  ijY*  +(i/,—  y,)"2 


V 
8, j  =  arctau 


.'/.  -  '/> 


Xi-  x, 


1.  Flow  Tangency  Condition 

The  flow  tangency  condition  in  the  case  of  no  blowing  or  suction,  is  satisfied 
at  the  middle  of  each  panel  sometimes  called  the  no  penetration  condition.  To  obtain 
this  condition,  differentiate  equation  3.10  with  respect  to  the  normal  vector  for  each 
panel. 


d-cixx,Ui) 
dn, 


=  H.7i,-  =  0 


v  v 

E  AAijffj  +  7  £  BBlJ  =  sin(fl,  -  a) 


.11! 


f    t  0 


12] 


where. 


,J       2zJ,  On,      v 

I    f  d              f  (.7.—  -V>)     ., 

dp     _ /         arctau    ; ««> 

""          2- J,  dm              [[*i-*j)l 
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AAij  ■  •    The  influence  coefficient   or  normal  velocity  (normalized  with   VQ)  at  the 
midpoint  of  the  i-th  panel  due  to  a  source  distribution  at  the  j-th  panel. 

BBij  •  •  •  "Fhe  influence  coefficient  or  normal  velocity  (normalized  with   VQ)  at  the 
midpoint  of  the  i-th  panel  due  to  a  vortex  distribution  at  thej-r/z  panel. 

xt  ,  //i  ■  .    Coordinates  of  panel  midpoint  of  i-th  panel 

x  .  ,  ,j  .  .    Coordinates  of  panel  midpoint  of  j-th  panel 

(j length  of  j-th  panel 

2.  Concept  of  the  Influence  Coefficient 

Eauation  3.12  can  be  evaluated  bv  intearation.  To  obtain  this,  we  have  to 
relate  the  coordinate  of  a  point  in  the  j-'th  panel  with  the  known  coordinates  of  the 
boundary  points  and  panel  angles.   From  the  geometry  in  Figure  3.2,  one  obtains 


2£i  =  c«J,  (-"I 


drii 
d'Ji 


212..  da  J,  (3-"l 


but,  cos#.  =  -  sin  9.  and  sin  J3.  =    cos  9.  ,  therefore  equation  3.13  and  3.14  become 

P-  =  -sm9i  (3.15) 


dtii 


|^  =    COS*  (3.16) 


On  the  other  hand. 


*j  =  Xj  +  tjCt*8j  (-I 


(3.1S! 


And  the  influence  coefficients  can  be  obtained  in  terms  of 
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^,j=h\cF-DG\ 
BB.j^^-DF-CG] 


3.19) 
3.20) 


where, 


A  =  -  ( x, -  X} )  cos  0j  -  ( .(/,  -  Yj )  sin  Bj 

B  =  [x,-X})2  +  [U>-Yj)'i 
C  =  sinj  9t  -  0j ) 
D  =  cos(0<  -  Q}) 

E  =  \/B  -  A-  =  ( Xi  -  A*, )  siu  <9;  -  ( .y,  -I"))  cos  fl; 

/."•  +  2.4/ 
F  =  In  1 


G'  =  arc ran 


5 
£7* 


L5-.-L/J 


1 


wnen  ;  =  j 


AAa  =  -         aud         BB,i  =  0 


3.  Computation  of  Total  Disturbance  Velocity  V 
The  computation  of  the  total  disturbance  velocity  in  x  and  y  direction  can  be  obtained 
bv  differentiating  in  the  x  and  v  direction. 


V  =  \'J  +  Vyj  =  grcidf 


(3.21) 
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.-L-ify     .    The  horizontal  velocity  component  at  midpoint  of  i-th  panel  due  to  a  unit 
source  distribution  at  the  j-th  panel. 

AAy      .   The  vertical  velocity  component  at  midpoint  of  i-th  panel  due  to  a  unit 
source  distribution  at  the  j-ih  panel. 

BBfj  .  .   The  horizontal  velocity  component  at  midpoint  of  i-th  panel  due  to  a  unit 
vortex  distribution  at  the  j-th  panel. 


SB* 


The  vertical  velocity  component  at  midpoint  of  i-th  panel  due  to  a  unit 

vortex  distribution  at  the  j-th  panel. 
These    influence    coefficients    can    be    obtained    by    integrating    equation    3.12    and 
introducing  equations  3.15  thru  3.18  into  equation  3.12 


AAfj=  ±£-lc<*ejF  +  sinfy3| 


(3.24) 


BBfj=  — [  -i-sin^F-cos^G'l 


(3.25) 


<U!j  =  BBr 


(3.26) 


BBlj.-AAfj 


(3.27) 
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when  i  =  j 


AAft  =  -  -  sin  9i 

3B?,=  |cos* 

AA»t  =  BBf, 
BB*t  =  -A4J 


4.  Kutta  Condition 

As  in  all  problems  concerning  airfoils  in  inviscid  flows,  an  auxiliary  condition 
needs  to  be  invoked  to  ensure  that  a  unique  solution  is  obtained.  This  condition, 
known  as  Kutta  Condition,  relate  to  assumptions  about  the  How  characteristics  at.  or  at 
least  in  the  neighbourhood  of,  the  trailing  edge.  When  the  surface  velocities  are  made 
equal  at  the  midpoints  of  the  trailing  edge  elements  then,  by  the  Bernoulli  Equation,  the 
pressures  at  these  points  are  also  equal,  so  the  Kutta  condition  can  be  represented  as 
the  condition  of  zero  loading  in  the  region  of  the  trailing  edge'  [Ref.  1]  which  is 
physically  realistic.    Therefore  the  Kutta  condition  can  be  summerized  as  follows: 

1.  For  a  given  angle  of  attack,  the  value  of  circulation  T  around  the  airfoil  is  such 
that  the  flow  leaves  the  trailing  edge  smoothly. 

2.  If  the  trailing  edge  angle  is  finite,  then  the  trailing  edge  is  a  stagnation  point. 

3.  If  the  trailing  edge  angle  is  cusped.  then  the  velocities  leaving  the  top  and  the 
bottom  surfaces  at  the  trailing  edge  are  finite  and,  equal  in  magnitude  and 
direction. 

These  are  the  basic  principles  of  the  Kutta  condition.  In  the  actual  computation  we  are 
using  in  the  code,  there  is  no  restriction  whether  the  trailing  edge  angle  is  cusped  or 
not  as  mentioned  in  item  3,  but. rather,  two  tangential  velocities  in  the  last  control 
points  assume  have  the  same  magnitude  but  in  opposite  direction. 

In  the  numerical  solutions,  the  actual  trailing  edge  is  not  a  stagnation  point. 
Furthermore,  it  is  found  that  the  velocities  at  the  midpoint  of  the  trailing  edge  elements 
differ  significantly  from  stagnation  values;  they  are  more  likely  to  be  closer  to  the  free 
stream  velocities. 

5.  Determination  of  the  vortex  strength 

The  Kutta  condition  is  used  to  determine  the  vortex  strength  y.  From 
equation  3-12  ,  set  y  =  1.0  and  obtain  the  expression  of  c  in  terms  of  the  influence 
coefficient 
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E    AAijfj  ml 


X 


E   BBij 

L  ;=i 


siu(#,  -  or) 


(3.28) 


This  equation  can~be  solved  by  using  Gaussian  elimination.  The  result  can  be 
substituted  into  equation  3.19  and  3.20  .  By  letting  /=  /  and  /=  A*  for  the  trailing  edge 
panel,  we  can  compute  horizontal  and  vertical  velocities  for  i-ih  and  ;V-//i  panel  as  I'  ,, 
I'  i,  r  y  and  K.y  Introducing  these  values  into  the  Kutta  condition,  produces  a 
quadratic  equation  in  y.  This  method  is  not  the  best  way  to  satisfy  the  Kutta  condition 
in  steady  How  calculation,  but  by  using  this  method  the  code  can  be  extended  to  the 
unsteadv  case. 
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<r  and  Pj  are  the  result  obtained  from  equation  3.2S.  From  equation  3.29.  there  are 
two  values  of  vortex  strength'  y,  but  only  one  value  is  used  in  the  further  calculation. 
These  values  relate"  to  the  angle  of  attack.  If  the  angle  of  attack  is  positive,  use  the 
negative  value  of  y.  If  angle  of  attack  is  negative  use  the  positive  value. 

6.  Determination  of  the  source  strength 

Once  vortex  strength  has  been  determined  using  equation  3.12,  the  source 
strength  <r  can  be  obtained  from 

W=  -K1  +  [J  A  (3-30) 

7.  Calculation  of  'on  body'  velocities 

The  velocity  at  midpoint  of  the  i-th  panel  can  be  obtained  by  a  spatial 
derivative  of  the  velocity  potential  in  the  tangential  direction. 

V;  =  i-[?(.r,,if.-)]  (3-D 

v  -v 

Vi=    E   .-Lri;<7;  +  7  E   B£7  +  coa(0,-a)  (3-32) 

where      ATij  =  -^[-CG-  \dF\  (3.33) 


BTij=U-DG+\cF\  (3-34) 


when  i  =  j      ATu  =  0.0  (3.35) 


3T"  ~  o  (3.36) 


AT, j  .is  the  influence  coefficient  of  tangential  direction  at  the  midpoint  of  the  i-th 
panel  due  to  a  unit  source  at  j-th  panel,  and  BT,j  is  the  influence  coefficient  of 
tangential  direction  at  the  midpoint  of  the  i-th  panel  due  to  a  unit  vortex  3.1  j-th  panel. 
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a.    Calculation  of  Pressure  Coefficients 

Once  the  velocity  on  the  body  has  been  calculated,  we  can  easily  obtain  the 
pressure  coefficientrby  the  Bernoulli  equation 


Pi  ~  PO  ,        t-2  n  ,7x 

—  =    1  -  V  ,  (j.J/) 


&? 


where  V.  is  the  dimensionless  velocity  (normalized  by  VQ)  at  the  midpoint  of  the  i-th 
panel.  Furthermore,  as  soon  as  the  pressure  coefficients  are  known  other  coefficients 
such  as  Lift  Coefficient  CL),  Drag  Coefficient(CD)  and  Moment  Coefficient  about  the 
leading  edge  (CM)  can  be  obtained  by  using  pressure  integration  along  the  body 
contour  which  is  approximated  by  a  closed  polygon. 


.v 

c*=  E  cp,.[(i-t+l-rt)i 


(3.38) 


i=i 


Cy=   E    -CPi.[(.Y(+1-A-,)I  (3-39) 


1  =  1 


.V 

CM=   E  CV[(-Y,-+i-  Xi)xi  +(I"i+i  -  YiM  (3.40) 

i=i 


CD  =  Cz  cos  a  +  Cv  sin  a  (3.41) 


CL  -  Cy  cos  a  -  C;  siii  a  f*  42) 


C.       LINEARLY  VARYING  VORTEX  DISTRIBUTION 

The  following  method  is  only  one  variation  of  the  use  of  the  panel  method 
[Ref.  7]  ;  it  involves  representation  of  the  airfoil  by  a  closed  polygon  of  the  vortex 
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panels.  The  vortex-panel  method  introduced  here  has  the  feature  that  the  circulation 
density  on  each  panel  varies  linearly  from  one  corner  to  the  other  and  is  continuous 
across  the  corner  as  indicated  in  Figure  3.4.  The  airfoil  and  wing  problems  can  be 
solved  by  means  of  a  vortex-panel  distribution  alone,  but  calculation  of  fuselage  and 
nacelle  characteristics  and  their  interference  flows  dictate  the  use  of  source  and. 
possibly,  double:  as  well  as  vortex-panels.  For  steady  flows,  it  has  been  found  that  the 
use  of  piecewise  constant,  discontinuous  distributions  of  voracity  can  lead  to 
inaccuracies  such  as  oscillating  values  of  the  voracity  en  successive  panels.  The  use  of 
a  linearly  varying,  continuous  distribution  of  vorticity  eliminates  this  problem. 

In  the  presence  of  uniform  flow  VQ  at  an  angle  of  attack  a  and  m  vortex  panels. 
the  velocity  potential  at  the  /-^'control  point  (xj.-.)  is  defined  as 


•i.Ui)  =  t0(-r«cosof  +  #sina)-     L     /     — —  arcrau    ; j  <i.9,  (  3.-U] 

;  =  i     Jj         -~  L  \X*  ~  XJI  J 

where  -   sj  =  7y  +  ( ~j+i  —  1j )  J-  '  0  (3.45] 


is  the  vortex  distribution  which  is  linear  along  the  panel  and  continuous  across  the 

boundary  points. 

The  panels.  N  in  number,  are  assumed  planar  and  are  named  in  the  clockwise  direction, 
starting  from  the  trailing  edge,  and  boundary  points  selected  on  the  surface  of  the 
airfoil,  are  the  intersections  of  continuous  vortexpanels.  The  (iV+7)  values  of '{■  at 
boundary  points  are  the  unknowns  to  be  determined  numerically.  The  condition  that 
the  airfoil  be  a  streamline  is  met  approximately  at  control  points.  The  boundary 
condition  requires  that  the  velocity  in  the  direction  outward  normal  vector  n.  be 
vanishing  at  the  i-th  control  point,  such  that 


r—  r[ri.Vi)  =  *        .  fori=  1.2.3....N  (3-46) 

0nt 


(A*,.  J-,)        Pan, 


nei  2 


Figure  3.4     Replacement  of  an  airfoil  bv  vortex  panels 
of  hncarlv  van/me  vortex  strength. 
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Carrying  out  the  calculation  of  equation  3.44  and  normalizing  the  vortex  strengths  by  2 
K  Vq  we  have 


E 


;=i   Jo 


dn. 


arctau 


Ijti  ~  fjj) 


rf«/  =  -sin(0t-  a)       (3-47) 


and  from  differentiation  and  intesration  we  set 


T  [CXitJ  f'j  +  (7*Vi#i  7';+i  I  =  sinM  -  or) 


j=i 


(3.48) 


where,   y'    =    y  -^/j  is   a   dimensionless   vortex   strength.     The   coefficients   in  the 
parentheses  are 

CXi    =  \dF  +  GG-GN*u 


where      .-L  =-(/,-  -  A"; )  cos  0,  -  ( m  -  I", )  sin  Bj 
B=  (x,—  Xj)2  +  [ui-  Yj)" 

C  =    WBL(&i-Bj) 

D  =  cosj  0t  -  9j ) 

E  =  (x,-  -  A*;)  sin  *,•-(#-  I*;)  cos fy 


G  =  arc  ran 
0 


£/ 


.  5  -  .4/j 
(*  -  Xj)  cos(4  -  2*y)  -  (z/.  -  Ii)sin(0i  -  20j 


The  expressions  in  the  parentheses  on  the  left  side  of  equation  3.48  represent  the 
normal  velocity  at  the  i-th  control  point  induced  by  the  linear  distribution  of  vortices 
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on  they'-:/?  panel,  called  influence  coefficient  of  normal  velocity  (except  for  i  =  m+7). 
When  i  =  j,  the  coefficients  have  simplified  values 

CX.iu  =  i 

which  describe  the  self-induced  normal  velocity  at  the  i-th  control  point. 
I.  Kutta  Condition 

In  order  to  apply  the  Kutta  condition  in  the  theoretical  analysis,  we  need  to  be 
more  precise  about  the  nature  of  the  flow  at  the  trailing  edge.  It  was  mentioned  before 
that  the  trailing  edge  can  have  a  finite  angle  or  it  can  be  cusped.  The  statement  of  the 
Kutta  condition  in  terms  of  the  vortex  sheet  is  y  (7"£)  =  V  -  V.  where  Vy  and  V j  are 
the  tangential  velocities  at  the  upper  and  lower  side  of  the  trailing  edge.  However,  for 
the  finite-angle  trailing  edge,  the  velocities  at  the  upper  and  lower  surface  have  same 
magnitude.  Thus  the  Kutta  condition  becomes  y  (TE)  =  0  , 

",[  +  Yjr+i  =  0  (3.49) 

For  computer  programming,  equation  3.48  can  be  arranged  such 'that  the- equation  is 
easv  to  formulate  in  matrix  form  •    ■ 


.v 

E   AXiM  =  RHSf        ,  i=1.2.3...,  ft+1  (3.50) 


After  combining  equation  3.49  and  3.50  they  are  sufficient  to  solve  for  the  (iV+/) 
unknown  */'•  values,  in  which  the  influence  coefficients  can  be  classified  as 

for  i  <  .V  +  l:      AXu  =  CXUl 

AX, ,  =  CXi     +  CXn  .  j  =  2.3 N 

j  i .  j  -,  j  j 

j  r  si  \m 

RHSi  =  mi(0;  -  a) 
for  i  =  .V  +  1:      AX a  -  .-LVi.v+i  =  1 

AXij  =  0  .  j  =  2.3 M 

RHSt  =  0 
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From  the  above  equation,  it  can  be  expressed  in  matrix  form 


AX-n 


A.\ 


.vi 


.■LVi.v     .4iVi.v+i 

.-L\.ViV+i 

0  1 


7i 
7i 


(    sinf^i  -  a)  } 

s'm(0-2  -  a) 


7a- 

7iV+l 


^  (3.51) 


sin((9.v  -  a] 
0 


To  obtain  the  unknown  */'.  values,  the  above  expression  can  be  solved  by  using  the 
Gaussian  Elimination  method  or  any  other  linear  equation  algorithm. 
2.  Calculation  of  the  'on  body'  velocities 

The  unknown  vortex  strength  having  been  determined,  we  can  proceed  to 
compute  the  velocities  and  pressures  at  every  control  point.  Recall  the  no  penetration 
condition  at  the  control  point.  Hence  the  disturbance  velocity  induced  at  control  points 
is  only  the  tangential  velocity.  The  velocity  can  be  obtained  by  differentiating  equation 
3.44  in  the  direction  of  the  tangential  vector  on  the  i-th  panel,  hence 


0<i 


3.52] 


r,  =cos(#,-a)+   E   [CT^  +  CT^'j+l]        .1=1.2.3 N       (3.53) 


in  which  the  coetlicients  in  the  parentheses  are  defined  as 


CTi     =  -CF-  DG-CT2ii 


1PF 


a 


CT-2 ...  =  C+  -  —  +  IAD-CE  \r 


2. j 


2   h 


lh 


when  i  =  j      CTUi  =  CTUi  =  -ff 


The  expression  in  the  parentheses  following  the  summation  symbol  has  the  physical 
meaning  of  the  tangential  velocity  at  the  i-th  control  point  induced  by  the  vortices 
distributed  in  the  j-th  panel.  Equation  3.53  can  be  simplified  further,  to  facilitate 
computer  programming, 
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.V 


Vt,  =  co<$(#,  -  a)  +   E   ^r,,7'      i  =  1.2.3 X 

and  AT,  j  is  defined  as      ATn  =  CTtl 

ATtJ  =  CTit}  +  CT2ij  j=  2.3, 

-^^"i.v-ri  =  cri,.jV 

in  convenient  matrix,  the  tansential  velocitv  can  be  written  as 
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3.55) 


3.  Calculation  of  Cp.CL.CD  and  CM 

After  the  velocities  at  control  points  have  been  determined,  one  can  obtain  the 
pressure  coefficients  at  the  i-th  control  point  using  the  Bernoulli  equation  3.-37  . 
Similarly,  other  coefficients  can  be  obtained  easily  by  pressure  integration  along  the 
bodv  contour. 


D.       DISCUSSION  OF  THE  PROGRAM  PANEL 

The  Panel  program  used  in  this  analysis  consists  of  two  FORTRAN  programs 
which  involve  the  implementation  of  source  and  vortex  distributions,  and  the 
implementation  of  linear  vortex  distributions.  Only  the  first  program  is  included  in  this 
thesis  (see  Appendix  A)  as  a  sample  program.  The  panel  methods  used  here  are 
basically  the  same  as  the  Hess  and  Smith  method  [Ref.  4]  where  the  source  is  located 
on  the  mid-point  of  each  panel,  constant  along  the  panel  but  different  from  panel  to 
panel  .  On  the  other  hand,  the  vortex  is  considered  constant  for  all  the  panels.  To 
satisfy  the  Kutta  condition  at  the  trailing  edge,  the  velocity  components  are  computed 
by  differentiating  the  velocity  potential  in  the  horizontal  and  vertical  direction.  And 
introducing  these  values  into  the  Kutta  condition,  we  can  solve  for  the  unknown  vortex 
strength  and  subsequently  for  the  source  strength. 
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START 


J 


'Read  input  file  code  "5" 
Number  of  nodes.  AOA, 
Airfoil  coordinates 


Compute  the  s 
for  each  panei 


ope. length. 


jJl 


Commute  the  influence 
Coefficients  for  normal, 
tansential.  x  and  y  direction. 
Setup  the  matrix  equations 


Set  y=  1.  solve  the  svstem  of 

equation  bv  usins  Gaussian  s 
Elimination  '.vithTtwo  RHS 


Meet  Kutta  Condition,  solve 
Vortex  Strength  r/)  and 
Source  Strength  ('en 


Compute  and  print  out  the 
Pressure  coefficients. 


_k_ 


Pressure  Intesration  to 
compute  CD.TTL  and  CM 


\/_ 


STOP 


Fisure  3.5     Source  and  Vortex  Panel  Method. 
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In  the  second  method,  only  vortices  are  used  as  singularity  distributions.  The  only 
difference  here  is.  that  the  vortices  are  distributed  linearly  along  the  panel  and 
continuous  from  one  panel  to  the  other. 

1.  Input  Data 

The  input  data  for  the  program  PANEL  must  be  arranged  in  the  following 
order: 

1.  Header  card.  The  header  card  consists  of  the  input  for  the  number  of 
coordinate  points  (column  1-10,  integer)  and  input  for  angle  of  attack  (column 
11-20  real) 

2.  Coordinates  of  the  airfoil.  The  arrangement  for  the  body  coordinates  must  be 
inputed  in  the  following  sequence:  start  from  the  trailing  edge,  progress  on  the 
lower  surface  to  the  leading  edge,  return  through  the  upper  surface  and  finish  at 
the  trailing  edge,  so  that  the  trailing  edge  coordinate  will  be  accounted  twice. 

2.  Program  output 

The  output  from  this  program  (see  Appendix  B)  can  be  arranged  as  follows: 

TABLE  1 
DESCRIPTION  OF  THE  PROGRAM  OUTPUT 

1.  PNL(I)  :  is  the  panel  number 

2.  X(I)  :  x  coordinate  of  control  points  (mid-points) 

3.  Y(I)  :  y  coordinate  of  control  points  (mid-points) 

4.  VEL(I)  :  is  the  dimensionless  velocity  for  each  control  point" 

(total  velocity  divided  by  free  stream  velocity) 

5.  GAMMA(I)  :  vortex  strength(must  be  the  same  for  the  each  panels) 

6.  SIGMA(I)      :  source  strength  for  each  control  point 

7.  CP(I)  :  pressure  coefficients  for  each  control  point 

8.  CL.CD.CM    :  coefficient  of  lift,  drag  and  moment  respectively 


E.       DISCUSSION  OF  THE  RESULTS 

Figure  3.6  shows  some  results  from  the  present  computation  for  various  angles  of 
attack.  In  general,  both  methods  give  very  good  agreement  in  pressure  distribution, 
although  there  are  slight  differences  in  the  region  of  minimum  pressure  and  at  the 
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trailing  edge.  In  the  leading  edge  neighborhood  ,  the  source  strength  may  give  a 
significant  effect  in  the  region  where  the  velocity  changes  rapidly.  On  the  other  hand, 
the  different  implementation  of  the  Kutta  conditions  at  the  trailing  edge  also  give  slight 
differences  in  the  pressure  distribution.  Recall  that  the  first  method  using  both 
velocities  at  the  first  and  last  panels  are  the  same  in  magnitude  but  in  the  opposite 
direction,  while  the  second  method  is  using  the  condition  of  zero  vortex  strength  at  the 
trailing  edge. 

From  the  above  conditions,  the  conclusions  can  be  obtained  from  the  present 
computations.  Both  methods,  can  be  used  to  predict  forces  and  moments  around  the 
airfoil  as  long  as  the  fluid  assumed  remains  inviscid  and  is  not  changing  with  time. 
Experience  dictated  that  the  second  method  gave  less  execution  time  than  the  first 
method  might  be  true  because  in  the  second  method  solve  only  one  unknown  (y) 
instead  of  two  (y  and  <x)  for  the  first  method.  On  the  other  hand,  the  second  method  is 
more  complicated  in  the  numerical  'formulation  than  the  first  one,  because  of  the 
varvina  vortex  strenath  throughout  the  airfoil  contour. 
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IV.  VISCOUS  FLOW  METHOD 

In  fluids,  momentum  is  transferred  by  internal  stresses,  namely  the  hydrostatic 
pressure  and  the  viscous  stresses.  When  fluids  are  affected  only  by  pressure  and  not  by 
viscous  stresses,  their  behavior  is  relatively  easy  to  predict  by  standard  inviscid  flow 
methods  as  described  in  chapter  3. 

A  variation  of  velocity  in  the  direction  normal  to  the  direction  of  the  velocity 
itself  is  called  a  shear.  Especially  in  high  Reynolds  Number  flows  this  shear  layer  is 
very  thin. 

The  most  common  type  of  a  shear  layer  is  the  boundary  layer  between  a  stream  and  a 
solid  surface.  On  the  solid  surface,  the  fluid  velocity  is  reduced  to  zero  (no  slip 
condition),  but  there  is  no  direct  constraint  on  the  velocity  gradient  at  the  surface.  At 
the  outer  edge,  the  velocity  tends  asymptotically  to  the  free  stream  values. 

As  mentioned  before,  in  1904  L.  Prandtl  came  up  with  his  well  known  theory  of 
boundary  layers.  Prandtl's  hypothesis  divides  the  flowfield  past  a  body  into  two 
separate  regions,  namely: 

1.  The  region  very  close  to  the  body  where  viscous  effects  are  important. 

2.  The  remaining  region  where  inertia  terms  are  more  dominant  than  viscous 
terms,  so  that  this  region  can  be  treated  as  inviscid  flow. 

These  assumptions  allow  us  to  deal  with  the  parabolic  boundary  layer  equations, 
instead  of  the  elliptic  Navier-Stokes  equations.  The  prior  experience  indicates  that 
parabolic  equations  can  be  solved  very  rapidly  and  efficiently.  Numerically,  the  change 
of  characteristics  means  a  change  from  a  field  procedure  to  a  marching  method,  which 
integrates  the  boundary  layer  equation  for  given  initial  conditions  step  by  step, 
proceeding  in  the  downstream  direction.  Depending  on  the  boundary  conditions, 
boundary  layer  methods  fall  into  three  types: 

1.  The  direct  boundary  layer  method.  This  method  employs  the  'no  slip 
condition','  requiring  zero  normal  and  zero  tangential  velocity  at  the  surface, 
and  a  prescription  of  the  external  velocity  at  the  edge  of  the  boundary  layer. 

2.  The  inverse  boundary  layer  method.  The  prescription  of  wall  shear  or 
displacement  thickness  replaces  the  above  edge  boundary  condition. 
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3.     The    interactive    boundary    layer    method.      The    edge    boundary    condition 
prescribes  a  combination  of  displacement  thickness  and  external  velocity. 
Further  discussion  "Will  focus  on  the  first  and  the  third  method,  since  the  computer 
code  uses  those. 

A.       DIRECT  BOUNDARY  LAYER  METHOD 

Direct  methods  are  used  in  the  region  where  the  viscous  effects  remain  small.  The 
current  code  applies  this  method  in  the  region  near  the  leading  td2s.  Direct  methods 
allow  the  generation  of  initial  conditions  at  the  stagnation  point  and  the  efficient 
integration  around  the  leading  edge.  The  numerical  approach  features  a  finite 
difference  method,  which  recasts  the  continuity  and  momentum  equation  as  a  system 
of  linear  algebraic  equations  [Ref.  13.]  To  begin  with,  we  consider  2-D.  steady  Hows  of 
incompressible  fluids,  described  in  a  curvilinear  coordinate  system  with  \  directing 
along  the  airfoil  surface  and  y  perpendicular  to  the  airfoil  surface.  The  velocity 
components  u  and  v  shall  be  determined  such  that  they  satisfy  anywhere  in  the  flow 
field  the  continuity  equation  (4.1)  and  the  momentum  equation  (4.2) 

dn       dr 


dn  dn  dn.  d  ,    On 

dx         dy         ''  dx         dy     dy 


11  7T  +  {'  —  =  «*  — ^i'—  (6— )  (4.2) 


where  3  boundary  conditions  are  needed  at  the  boundaries  of  the  flowfield, 

U  =  0:         u[x.O)  =  0.  c(x.O)  =  0 

y  =  ye  :       u{x.i/e)  =  at[x) 

with  b  =  1  +  V  ,  V.  These  equations  are  referred  to  as  boundary  layer  or  thin  shear 
layer  equations.  To  solve  these  equations,  it  is  convenient  to  introduce  a  stream 
function  (  u  =  c\\i  dy  and  v  =  -  c\\l  dx  ).  which  reduces  the  number  of  dependent 
variables.  Since  the  stream  function  automatically  satisfies  the  continuity  equation, 
only  the  momentum  equation  is  left 

dv     d'u        dc   d'2i.'  dn.  d       d'v 

—  — t  =  n. — -  +  v—{b 1  (4.3) 

dy   OxO'j       ox    dy-  '  dx  dy     dy 
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This  equation  is  subjected  to  the  Falkner-Skan  transformation,  which  scales  the  normal 
coordinate  y  and  the  stream  function  \\t  with  reference  to  the  external  velocity. 

t'r- y  " 

ux 
/(-'/)=  -=tt-(x.//) 

With  these  transformations,  the  momentum  equation  becomes 

f*n'+2~±/r+»[i-  (/')»!  =  ,(/'^-/»^)  ,4.4, 

\     ox  dx* 

and  the  boundary  conditions  are 

-7  =  0:       /%-.())    =Q.      /(...o)  -0 
'/  =  '/*:       /'(x.'/«)  =  1 

where  m  is  defined  as  a  dimensionless  pressure  gradient  parameter  !  m  =  (x  u  )(du  dx) 
}.  and  prime  denotes  differentiation  with  respect  to  r\.  This  is  a  third  order  partial 
differential  equation.  The  solutions  of  this  PDE  are  called  non-similar  flows,  because 
they  depend  on  both  x  and  r\.  In  contrast,  if  the  right  hand -side  of  equation  4.4 
vanishes,  and  therefore  the  solutions  depend  on  r|  only,  they  are  known  as  similar 
Hows. 

1.  The  Box  Method 

One  of  the  most  flexible  methods  in  solving  non-linear  differential  equations  is 
the  box  method  developed  by  Keller  [Ref.  13]  in  1970.  The  basic  steps  of  the  box 
method  are  the  conversion  of  the  governing  equations  into  a  first  order  system,  the 
discretization  of  the  differential  equation  by  using  central  differences  and  two  point 
averages,  the  linearization,  and  the  solution  of  the  resulting  algebraic  system.  The 
introduction  of  two  additional  dependent  variables  U  and  V  converts  the  third  order 
momentum  equation  (4-4)  into  a  first  order  system 

f'  =  U  (4.5) 


U'=V  (4.6l 
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(»■)'+ £?V+«U-I7»|  ~*{U%L-Y9-L\  (4.7) 

and  the  boundary  conditions  become 

7  =  0:       ^(;.0)   =0.     /(;,(.])  =  o 

Instead  of  dealing  with  continuous  functions  /,  V,  and   V,  we  use  a  set  of  discrete 

values  of  these  flow  properties.    Let  the  solution  domain  0  <  i  <  //.    0  <  /;  <  ;/j     be 

covered  by  a  rectangular  mesh  (Fieure  4.1) 

x\  =  0.         X,-  =  Xi-i  -  ^      with     2  <  i  <  I 

//x  =  0.  i] j  =  r/;_i  -i-  /j;      with     2  <  j  <  J     and     qj  -  qe 

We  approximate  all  quantities  whether  it  is  a  function  or  its  derivative  or  a  parameter 
like  b'.  in  terms  o[  nodal  values  and  coordinate  of  the  network.  The  stream  function 
and  its  first  and  second  derivatives  with  respect  to  r|  are  abbreviated  at  the  nodes  of 
network  by 

{An.nj).  vUi-nj).  V[zirnj)}  =  If),  U).  v;} 

The  solution  of  the  parabolic  boundary  layer  equation  at  a  certain  streamline  position.0 

say  X;  depends  solely  on  the  solution  of  upstream  positions  say  x-  j  .  x-  >  .  .  .    while 

no  downstream  influence  has  to  be  considered.    The  overall  solution  can  be  obtained 

step   by   step  with   the  calculation   propagating   from  the   stagnation   point   into   the 

downstream   direction.     The   advantage    of  using    first    order   equations    and    central 

differences  is   that   we   can   reduce   the   domain   of  dependence   from 'all   upstream  x- 

stations  to  the  immediate  preceding  one.  hence  one  step  of  the  solution  procedure 

writes  the  governing  equations  for  a  column  of  net  rectangles  (boxes)  residing  in  the 

subdomain 

■Pr-i  <  '  <  J?i     and     0  <  rf  <  rjj 

and  solves  subsequently  for  the  nodal  values  of  the  downstream  face  of  the  rectangular 
shaped  subdomain.  The  x-station  being  currently  solved  holds  therefore  the  superscript 
"/",  while  "/-/"  denotes  the  known  flow  properties  of  the  adjacent  upstream  location.  As 
indicated  by  the  term  central  differences  the  equations  are  satisfied  midway  between  the 
nodes. 

The  two  ODE  's  are  centered  about  ( /,-,  iiJ_l/2) 

And  the  PDE  is  centered  about  (/,-i/2- '/;-!/::) 
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Figure  4.1     Net  rectangle  for  finite  difference  approximation. 
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So  equations  4.5  .4.6  and  4.7  can  be  written  in  terms  of  finite  differences 


hj        -J^  +  ^-i) 

J7J  -  cr«_f      i     . 

— -=  in    -i-  r*    ) 


Ibvyr^-ibWZ?        m'-Va  +  1 


(/n;:$  +  m'-v* 


-  *i-l/2 


r/ 


—  1/2    U.<-U-2        C'-i/2  T-«-l/2    A-l/2  ~   *.- 


7-1/2 


A", 


7-1/2 


(4.S) 

(4.0) 


4.10] 


and  the  boundarv  conditions  take  the  form 


Ul  =  0, 

vl  =  i 


A'  =  o 


The  subdomain  under  consideration  consists  of  J- 1  net  rectangles,  the  flow  quantities 
of  each  being  related  by  a  momentum  equation.  The  equations  4.S  and  4.9  link  the 
dependent  variables  to  their  rj-derivatives. 
2.  Newton's  Method. 

Unfortunately,  the  unknowns  appear  in  nonlinear  combinations.  Therefore  we 
introduce  Newton  's  Method  to  solve  this  nonlinear  system.  The  solution  involves  an 
iterative  procedure,  in  which  the  variables  are  linearized  around  their  values  of  the 
preceding  iteration 


j  j      j  j 


y,.K 

J 


<.K-1 
I.K-L 
.K-l 


^ 


.i.i 


j  j  j 


for  k>  2.     U'A  =  C'!"1  +  tt 


+  «■;■« 


'    J  J 


where     5fj     <  / 
where     6U''K  <  U 


t.x-i 


I.K-i 

7  "7 

•  •  K      ^*    T  ' '  •  *—  1 


where     51" '    <.  V 


for  «  >  2. 

where  k  denotes  the  iteration  counter.  Substitution  of  these  values  into  equations  4.S. 
4.9  and  4.10,  and  dropping  the  quadratic  terms  in  (<!>/''*,  6U%'K,  5V*'K  )  lead  to  a 
linear  system  in  the  unknowns    i>f'K.  6U1,'*.  5\"'K 

7  7 
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"!*  -  «/m  -  t  ( «tf"  +  «fc )  =  /;.Tl  -  /r1  +  hj  vfi*      , 4.11 ; 


where  the  coefficients  in  the  bracket  are  defined  as 

;'.AC-i 


'  hj  2kt      l/^-V2       -0-1/2  J  4      /;-l 


»3 


V 


2A\     (    >-V2        ';-V2>    '        4       v; 


ffl'^1    T-l'./C-l 


l**Jj      -      0,         *   \  -1/2    "S-i/2)  ~  ~ ^ ^ 


,f,r.-(_iu.  +  Mi)tr; 


r^);   -- j j- +  -y-(/»  )>-i/2  -  ( — 3^; —  +  m  M' :  j-i/2' 

f '- t/2     /T-I.K'—  1     fi.K—l  T-«-l        fi.K—  1  ri-l        T-l'.K— l\l 

\     j-lj'l   JJ-U'i    "*■      J- V2  ;7-l/2         7;-l/2  S-l/2  'J 


\  h}  2        (/;  ';-i/2  +  (— ^ "'       )([Vi/2) 


:0 


3.  Keller's  Block  Elimination  Method 

Together  with  the  boundary  conditions  this  system  is  repeatedly  solved  until 

6ft'K.5U''K.5Y''K  ~are  small  enough  to  be  neglected.    Equation  4.11,  4.12.  and  4.13 

can  be  solved  by  Keller's  Block  Elimination  Method.    Block-tridiagonal  matrices  are 

composed  of  submatrices.  called  blocks,  of  which  only  those  residing  on  main  and  both 

adjacent  diagonals  have  non-zero  entries. 


•"4    i     i°t    1 


;!     \c\ 


r  /-»<.« 


Wljl  U/iJ  [ci-ii 


-V  i 


l'-2     J 


to 


{')■"}    =    W 


{r'jl 


(4.14) 


*'her?  the  biocks  are  3  x  ^-matrices  bein?  de£aed  as 


-Ay/2  0 


'3jr  (*5);-x 


•»1. 


"I  -/l;+l/2 


for  2  <  y  <  /  -  1 


[£?;•«]  = 


-1  -hj/2         0 

My*  M';K   (f:Y;K 

0  0  0 
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1  0  0 
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0     -1  -/ia/2 
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for  2  <  ;  <  J 


for  1  <  y  <  /  -  1 
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and  where  the  ri>ht  hand  sides  are  obtained  fro 


m 


/.K       fas  given  by 

lraJy    -1  .  for2<  j<  J 

J         { momentum  equation  ™  J  — 

( r: )  j-  =  er;~*  -  cr^  +  A/+l  i££  for  i  <  y  <  /  - 1 

Ml"* -0,  (rsJx'-Q.         (r3)j*  =  0 

The  unknowns  of  the  linear  equations  are  the  Newton  iterates  of  the  stream  function  ( 

/•**  ),  its  first  derivative  (  UJ>*  )  and  its  second  derivative  (  V-^  ).  This  method  is 
very  effective  and  it  consists  primarily  of  two  steps: 

1.  The  forward  step  eliminates  the  lower  diagonal  of  submatrices. 

2.  The  backward  step  solves  the  remaining  system  from  bottom  to  top. 

B.       INTERACTIVE  BOUNDARY  LAYER  METHOD. 

The  application  of  the  direct  method  is  restricted  to  regions  where  the  viscous 
effects  remain  small.  Integration  of  the  boundary  layer  equations  will  break  down  at 
the  point  of  zero  skin  friction.  To  avoid  this  break  down,  we  need  a  method  that  is 
able  to  integrate  the  boundary  layer  equations  through  the  point  of  zero  skin  friction. 
Further,  these  methods  are  required  to  account  for  strong  interaction  effects,  arising 
from  boundary  layer  separation  or  the  rapid  flow  acceleration  downstream  of  the 
trailing  edge,  both  of  which  cause  substantial  changes  in  the  external  velocity 
distribution. 

In  contrast  with  direct  and  inverse  methods,  the  interactive  method  treats  the 
external  velocity  and  displacement  thickness  as  unknown  quantities,  reflecting  the 
elliptic  character  of  the  outer  flows.  This  introduces  apparently  one  additional 
unknown  into  the  viscous  flow  problem,  whose  solution  can  be  obtained  by  using  two 
methods  : 

1.  The  eigenvalue  method,  or 

2.  The  Mechul  function  method. 

The  second  method  is  being  preferred  here,  since  the  first  method  involves  non-linear 
eigenvalue  problems.  The  edge  boundary  condition  of  the  direct  problem  is 
supplemented  by  the  so-called  interactive  boundary  condition,  which  relates  the 
unknown  external  velocity  with  its  inviscid  and  "displacement-perturbation-related" 
contributions.  Boundary  layer  equations  in  the  following  constitute  a  system  in  the 
unknown  functions  u(x.y),  v(xy)  and  ue{.x.y) 
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dn        d>' 

Ox         d<j~'lidx    ~U~d\}     ~j  (*-16J 


0--  (4.17) 


These  equations  consist  of  continuity  equation  (4.15)  ,  momentum  equation  in  x-dircction 
(4.16)  and  the  seemingly  unnecessary  momentum  equation  in  y-direction  (4.17).  where 
the  pressure  term  has  been  expressed  in  terms  of  the  external  velocity.  The  Mechul 
function  approach  assumes  that  the  external  velocity  be-  a  function  of  two  arguments, 
resulting  in  the  need  for  the  trivial  y-momentum  equation.  The  reason  for  considering 
u  (x1y)  rather  than  u  (x)  is  for  purely  numerical  reasons,  i.e.,  such  a  provision  allows  an 
easy  setup  of  the  finite  difference  equations  avoiding  the  eigenvalue  technique. 

The  govenng  equations  are  complemented  by  proper  boundary  conditions.  The 
velocity  components  u  and  v  are  required  to  satisfy  the  no-slip  conditions  at  the  surface 
of  the  airfoil  and  the  horizontal  component  must  merge  smoothly  into  the  outer  How 
at  the  boundary  layer  edge. 

i/=  Q:         «(.r.O)  =0.      ri-.O)  =  0 

>J=Uc-      «(*.&)  =  ««(*.J«),     ««(*•&)  =  »€/'(■*)  +  Z  I  m^u*S"^  JZT 

with  uAx)  denoting  the  inviscid  velocity  distribution  and  the  second  term,  called  Hilbert 
integral,  approximating  the  perturbation  velocity  due  to  viscous  effects.  The  interactive 
method  can  be  applied  to  attached  and  separated  How  regions,  while  direct  method 
cannot  do  so  because  of  their  breakdown  related  to  the  Goldstein  singularity,  nor 
inverse  methods  because  of  their  poor  convergence  rates.  Therefore  the  interactive 
methods  are  being  preferred  on  the  main  parts  of  the  airfoil.  Only  at  the  stagnation 
point  this  method  cannot  be  applied. 

The  steps  which  turn  the  partial  differential  equations  of  the  interactive  problem 
into  a  linear  system  of  algebraic  equations  resemble  those  of  the  direct  method,  so  that 
only  major  steps  and  their  results  will  be  repeated  here.  After  the  introduction  of  a 
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stream  function  {u  =  c\\f  cy  and  v  =  -c\\i  dx).  the  equations  undergo  a  transformation, 
which  scales  the  normal  coordinate  y,  stream  function  \\f.  and  the  external  velocity  u 
with  reference  to  the"  constant  velocity  uQ  and  the  local  streamwise  coordinate  x 

I «,) 

V   VX 
^U,)UX 

u.li.y) 

n  i*-n)  = 

where  uQ  is  taken  as  the  free  stream  velocity.  The  concept  of  constant  boundary  layer 

thickness,  attained  bv  Falkner-Skan  variables  with  u    as  reference  velocitv,  has  to  be 

e 

abandoned  because  the  external  velocity  is  unknown  in  the  interactive  calculations. 
Provided  that  the  integration  of  the  boundary  layer  equations  does  not  start  ;n  the 
immediate  neighborhood  of  the  stagnation  point,  the  growth  of  the  boundary  layer 
thickness  can  be  kept  limited.  In  terms  of  these  so  called  semi-transformed  coordinates 
the  boundary  layer  equations,  written  as  a  first  order  system  by  means  of  two 
additional  dependent  variables  U  and  I'.  and  the  boundary  conditions  take  the  form  . 

r  =  v  (4.is) 


U'  =  V 
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^"  =  0  (4.21) 


n  =  Q:        Ulx.O)   =  0.     /(/.0)  =  0 
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The  discretization  of  the  flow  field  follows  closely  the  above  outlined  procedure  of  the 
direct  method,  covering  the  generation  of  an  orthogonal  grid  and  the  introduction  of 
central  differences  and  two-point-averages.  The  overall  solution  proceeds  in  the 
downstream  direction,  accountinc  for  downstream  travelling  disturbances  only.    On  the 
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assumption  that  backflow  velocities  are  comparatively  small,  a  stable  integration  can 
be  carried  out  by  adopting  the  FLARE  approximation  (FlugQe-Lotz  and  Reyhncr).  The 
purpose  of  a  FLARE  approximation  is  to  permit  the  use  of  a  downstream-marching 
algorithm  in  regions  of  backflow.  This  is  accomplished  by  setting  the  streamwise 
convection  term  u  cu  dx  equal  to  zero  in  regions  of  backflow.  With 
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designating  an  " on- off  switch"  of  the  streamwise  convection  term,  the  finite  difference 
equations  of  the  interactive  boundary  layer  problem  become 
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Boundary  conditions  are  expressed  in  terms  of  nodal  values,  whereby  the  evaluation  of 
the  integral  occunng  in  the  interactive  boundary  condition  involves  an  approximation 
in  the  fashion  of  the  panel  method  approach,  leading  to 

U[  =  0.  /{  =  0 

where  £j  and  c^  denote  a  parameter  and  the  diagonal  element  of  the  interaction  matrix, 
resulting  from  a  discrete  approximation  to  the  Hilbert  integral.  Averaging  as  well  as 
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centering  is  supposed  to  obey  the  principle  that  the  number  of  generated  terms 
approaches  a  minimum.  This  entails  ordinary  difTerential  equations,  like  the  y- 
momentum  equation,  being  centered  about  the  middle  of  the  downstream  face  and 
partial  difTerential  equations,  like  the  x-momentum  equation,  being  centered  about  the 
middle  of  the  box. 

A    balance    of  unknowns,   which   occur   here   as    vectors   in   four   components. 
{/!•  C^,T»"'.  71"'} "confirms  the  principal  solvability  of  the  system.  The  J  quadruplets  of 
unknowns  match  with  a  total  of  4J  equations,  including  2(7-/)  auxiliary  relations.  (J- 1) 

x-momentum  and  (J- 1)  y-momentum  equations,  each  of  which  corresponding  to  one  of 
the  (J- 1)  net  rectangles,  and  4  boundary  conditions.  After  linearizing  this  system 
around  the  vaiues  of  the  preceding  iteration  (iteration  counter  "k-/"),  respectively 
around  the  solution  of  the  adjacent  upstream  x-station  in  case  of  the  first  iteration,  we 
arrive  at  a  linear  system  in  the  yen-ton  iterates  6  /'**.  5Ul'K.  6* ",'K.  ^n"'"* 
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supplemented    by    the   two    components    of  no-siip   condition,    edge   and    interactive 

boundary  condition 

yu"{*  =  0.      $f{'*  =  0 

6U'yK  -  6U'i'*  =  0 
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Terms  have  been  grouped  such  that  known  quantities  reside  on  the  right  hand  sides, 
while  unknown  quantities  appear  left  of  the  equal  sign.  The  abbreviated  coefficients  in 
the  momentum  equation  are  defined  by 
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and  the  right  hand  side  of  equation  (4.23) 
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Since  the  overall  procedure  involves  a  repetitive  linear  pattern  to  approach  the  solution 
of  the  nonlinear  system,  the  linear  equation  solver  has  to  be  as  fast  as  possible.  Fast 
algorithms  take  advantage  of  specific  matrix  structures,  one  of  which  is  the  block- 
tridiagonal  structure,  which  can  be  enforced  in  this  method  by  the  following 
arrangement  of  equations  and  boundary  conditions  : 
1.      First  set  of  equations  (index  y  =  /  ) 

1.  Wall  boundary  condition  prescribing  no  penetration 

2.  Wall  boundary  condition  prescribing  zero  tangential  velocity 
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3.  Second  auxiliary  relation  (linking  I"  and  V)  of  the  first  box 

4.  Trivial  y- momentum  equation  of  the  first  box 

2.  Intermediate  sets  of  equations  (index  2  <j<J-l) 

1.  First  auxiliary  relation  (linking/  and  U)  of  the  (J-I)-sr  box 

2.  Momentum  equation  of  the  (j-I)~st  box 

3.  Second  auxiliary  relation  (linking  I"  and  H  of  the  j-th  box 

4.  Trivial  y-momentum  equation  of  the  j-th  box 

3.  Last  set  of  equations  (indexy'=7) 

1.  First  auxiliary  relation  (linking/  and  U)  of  the  last  box 

2.  Momentum  equation  of  the  last  box 

3.  Interactive  boundary  condition 

4.  Edge  boundary  condition. 

Following  these  instructions  equations  and  boundary  conditions  can  be  arranged  in 
matrix-vector  form. 
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irhere  {#}•*}  =  {*/;•*,  «r;«,  «••;•",  ,in-;--}r   and        (r;:*}  =  {[rjf,^)';*,  (rzY;\ 
{i"-,)';*}7       are      f°ur  dimensional  vectors  of  the  unknown  Newton  iterates  and  the 
known  right  hand  sides,  respectively.  The  blocks  in  the  three  diagonals  of  the  above 
matrix  are  4x4  and  can  be  obtained  from 
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The  comoonents  of  the  risht  hand  side  vectors  can  be  calculated  from  the  following 


formulae 
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The  numerical  solution  of  the  above  system  can  again  be  achieved  by  Keller's 
block  elimination  method,  which  works  very  much  like  Gauss  s  algorithm,  but.  firstly, 
matrices  are  eliminated  instead  of  scalars.  and.  secondly,  quite  a  few  manipulations  can 
be  saved  because  of  sparse  occupation. 
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C.       INTERACTION  MODEL 

The  interaction  model  refers  to  the  coupling  of  the  boundary  layer  and  the 
external  inviscid  flow.  In  principle,  the  interaction  consists  of  thickening  the  effective 
airfoil  shape  by  viscous  displacement,  which  will  result  in  a  change  of  the  surface 
pressure.  Numerically  the  interaction  between  \iscous  and  inviscid  flow  regions  takes 
place  via  the  boundary  conditions  only,  which  is  the  specification  of  either  an 
impermeable  displacement  surface  or  a  nonzero  wall  transpiration  at  the  original 
surface  in  case  of  inviscid  flow,  respectively,  the  prescription  of  either  pressure,  or 
displacement  thickness,  or  a  linear  combination  of  both  in  case  of  viscous  flow.  If  the 
viscous  effects  on  pressure  remain  small,  then  the  interaction  is  called  weak.  However, 
situations,  where  viscous  disturbances  to  the  inviscid  flow  field  are  substantial,  demand 
the  application  of  strong  interaction.  In  general  .  interaction  models  can  be  classified 
into  four  types,  whereby  the  first  three  types  provide  loose  coupling  of  viscous  and 
inviscid  regions: 

1.  The  direct  interaction  method  combines  a  direct  inviscid  and  direct  viscous  flow 
solver.  This  classical  approach  achieves  a  solution  by  iterative  matching,  of 
boundary  layer  calculations  with  prescribed  pressure  and  inviscid  calculations 
with  prescribed  displacement  thickness.  The  alternate  treatment  of  pressure  and 
displacement  thickness  leads  to  a  local  breakdown  of  the  procedure  when  slight 
changes  in  pressure  entail  significant  changes  in  displacement  thickness  (see 
Figure  4-2a). 

2.  In  the  inverse  method  the  roles  of  the  displacement  thickness  and  pressure  are 
interchanged.  Hence  the  inviscid  flow  equations  determine  the  displacement 
thickness  distribution,  which  is  imposed  as  boundary  condition  on  the  viscous 
flow  calculation.  The  result  of  which  is  the  pressure,  which  closes  the  cycle  by 
being  input  to  the  inviscid  flow  computation.  The  hierarchical  manner  of 
solving  the  complete  flow  problem  excludes  this  model  just  as  the  previous  one 
from  handling  strongly  interacting  regions  (see  Figure  4-2b). 

3.  The  semi-inverse  interaction  method  is  composed  of  a  direct  inviscid  and  an 
inverse  visous  flow  solver.  Both  parts  of  the  scheme  process  the  same  input,  i.e. 
displacement  thickness,  and  generate  the  same  output,  i.e,  pressure.  After  each 
cycle  an  updated  displacement  thickness  is  relaxed  based  on  the  deviation  of  the 
two  pressure  distributions,  which  should  coincide  upon  convergence  (see  Figure 


65 


Airfoil 


Geometry 


v        INVISCID  (direct)       .£ 


*V) 


«eU) 


VISCOUS  (Direct) 


*"U) 


a)  DIRECT  METHOD 


VISCOUS  (inverse) 


n 


"  \x\ 


*M 


Airfoil 


Geometrv 


(b)  INVERSE  METHOD 


Airfoil  Geometrv 


INVISCID"  (direct) 


*•[*) 


»e/M 


1 


VISCOUS  (inverse) 
^ 


»«r(-r 


'*(/J 


RELAXATION 


*"(*) 


c)  SEMI-INVERSE 


Fieure  4.2     Direct,  Inverse  and  Semi  Inverse  method. 
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4-2c).  The  existence  of  a  definite  hierarchy  between  boundary  layer  and  the 
outer  flow  in  the  above  mentioned  methods  gives  rise  to  a  limited  rate  of 
feedback  between  both  regions.  Strong  interaction  models  incorporate  the  outer 
flow  somehow  in  the  boundary  layer  calculations,  for  example  by  the  following 
interaction  law: 
4.  Simultaneous  or  strong  interaction  methods  solve  the  boundary  layer  equations 
subject  to  an  interaction  law  as  outer  boundary  condition,  which  retrieves  the 
elliptic  character  of  the  outer  flow.  No  definite  assignment  of  displacement 
thickness  and  pressure  can  be  made  to  viscous  and  inviscid  region.  Rather,  both 
quantities  are  treated  as  unknowns,  related  by  the  interaction  law.  The 
procedure  emphasizes  simultaneous  solution  for  both  displacement  thickness 
and  pressure  (Figure  4. 2d). 

The  current  interaction  law  is  formulated  in  terms  of  displacement  thickness  5  and  the 
external  velocity  u  .  the  latter  being  related  to  pressure  by  Bernoulli's  equation.  The 
solution  of  the  boundary  layer  equations  in  an  iterative  fashion  makes  use  of  an  outer 
boundary  condition,  in  which  the  total  external  velocity  is  written  according  to 

ue(x)  =uellx)  +  ue6(x)  (4.31) 

where  u  Ax)  is  the  inviscid  external  velocity  and  us  t.x)  is  the  perturbation  due  to 
displacement  effects.  L'sing'  the  thin  airfoil  concept,  the  perturbation  velocity  can  be 
written  as  the  Hilbert  integral 

/    ^        l     f*m*M   ac       l    ft  i      ?•<     d* 


The  contribution  due  to  viscous  effects  can  be  evaluated  by  means  of  the  blowing 
velocity  concept.  Lighthill  proved  that  the  effect  of  boundary  layers  on  the  outer  flow 
can  be  represented  by  a  surface  distribution  of  sources.  The  source  strengths  must  be 
determined  such  that  the  virtual  displacement  surface  becomes  a  streamline  (see  Figure 

dx  uAx\ 
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Fiaure  4. 2d    Viscous-Inviscid  Interaction  Method. 
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where  vi.r.o  i  denotes  the  vertical  velocity  on  the  displacement  surface.  Since  it  is 
sufficient  to  approximate  the  correction  term  u.z  we  can  use  the  thin  airfoil 
approximation  as  follows: 

1.  The  uppef~and  lower  surfaces  of  the  airfoil  will  be  considered  as  flat  plates. 
This  implies  that  the  blowing  velocity  v(x,0)  equals  half  oC  the  local  source 
strength 

2.  The  displacement  thickness  is  assumed  to  be  so  small  that  the  horizontal 
velocity  components  of  the  inviscid  flow  do  not  vary  across  the  boundary  layer 

Oil)         '.  f6'dr    .  d 
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Jo    d'J  <k  (4.34) 


The  interaction  region  is  limited  to  the  finite  region  n,  <  z  <  /4on  either  the  upper  or 
lower  surface.  The  numerical  implementation  of  the  interaction  law  requires  some 
discrete  approximation  of  the  above  thin  airfoil  integral.  Adopting  the  panel  method 
approach,  which  concerns  here  a  piecewise  approximation  of  the  connnuous  blowing 
velocity  d{  ■c.S*)/dx  io  allow  piecewise  analytical  integration,  the  integral  can  be  written 


-  /      T7  M     7  =    I-  CikinJ  )  (4.35) 


wi:h[ct;.idenoting  a  matnx  of  interaction  coefficients  defining  the  relationship  between 
the  boundary  layer  thickness  and  the  external  velocity.  Recalling  that  the  boundary- 
layer  calculation  at  a  streamwise  position  involves  Newton's  iterations,  the  inviscid 
contribution  can  be  included  in  the  total  external  velocity  of  the  previous  Newton's 
iteration,  leading  to  a  generalized  version  of  the  interaction  law 

(u€)'-«»(iif)'-«-1+  E  Ctk{{u,y)^-{uj')k-*-1}  (4•:•6, 

Since  displacement  thickness  does  not  belong  to  the  dependent  variables.  (O'**.  the 
displacement  thickness  of  the  streamwise  position  currently  being  solved,  must  be 
expressed  in  terms  of  dependent  variables  to  make  allowance  for  its  unknown  status. 
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Ficure  4.4     Application  of  the  Direct  and  Interactive  Method. 
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With  (6*y-K    being  replaced  by      Uj  ~ty'*/( ue)'-K    and  after  separating  known  and 
unknown  terms,  the  interaction  law  takes  the  form 


[ue),'K(l-cuyj)  +  ciiVjK  =  gf  (4.3; 

with 

g<  -(«,)«.«-!  +   r  ^[ko*-*  -  (MY*'*"1]  -  c,,(«e^  -  or*"1 


Implementation  of  this  relation  necessitates  a  known  right  hand  side,  which  can  be 
evaluated  only  in  approximate  fashion,  because  the  term  (hj.S*)*"*  is  not  not  known 
yet  downstream  of  the  current  x-location.  To  ensure  interaction  also  over  this  region, 
these  terms  are  taken  from  the  previous  iteration  updated  by  some  relaxation  formula. 
With  dik  y/i/x/u0  and  gf  —  gf/u0  denoting  the  dimensionless  interaction 
coefficient  and  right  hand  side,  respectively,  the  actually  coded  interaction  law  can  be 
written  in  terms  of  semi-transformed  coordinates 


~K 


n-;-*(i-ct,///)  +  c„/;  =?r  (4.3s> 

This  equation  is  being  used  as  an  outer  boundary  condition  in  the  viscous  flow 
computation,  and  relates  the  unknown  external  velocity  with  the  unknown 
displacement  thickness  of  the  i-th  boundary  layer  station,  whereas  the  displacement 
thickness  has  been  expressed  by  the  streamfunction  and  external  velocity.  Because  of 
the  elliptic  character  of  the  outer  flow,  which  has  been  incorporated  in  the  boundary 
layer  via  the  interaction  law,  the  solution  requires  a  global  iteration,  consisting  of 
several  viscous  sweeps  to  be  performed  over  both  the  streamwise  upper  and  lower 
surface. 

Let  the  continuous  function  "external  velocity  times  displacement  thickness", 
denoted  by  D,  be  discretized  at  a  finite  number  of  streamwise  positions,  then  the  thin 
airfoil  integral  can  be  approximated  by  a  finite  series  of  weighted  "D's"  at  the  very 
locations 


lKdD     d£  L 

~rc r  =    2-  en,  Dk 


(4.39) 


The  weights  are  the  interaction  coefTicients  c-^  ,  with  the  first  subscript  indicating  the 
streamwise  position,  where  the  correction  term  u  s  is  to  be  evaluated,  and  the  second 
indicating  the  location,  whose  efTect  is  accounted  for.  With  a  properly  interpolated  D- 
function.  integration  can  be  carried  out  analytically  piece  by  piece.  Provided  the  point 
under  consideration  does  not  fall  within  the  limits  of  integration.  D  will  be  interpolated 
linearly.  A  piecewise  linear  function  can  be  built  up  by  overlapping  triangular 
distributions,  integration  over  which  yields  the  coefficients  whose  k*i ,  kti-l  ,  k*  i+  1 
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Cik  = 


(•■»**  i 


dD      di 


D, 


■i      dD 
with     — 

.di 


J?fc-Xfc-1 


Cik  -  - 


Iu! 


for    Xk-i  <  i  <  *k 
for    xk  <  i  <  Jfjfc+i 
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xk-xk-i 


x,-xk.-L 
Xi-Xk 


lii 


Xi-Xk 


Ik  +  l-  Xk         I  J*«  —  "Tfc+1   I. 


4.40) 


4.41' 


A  linearly  interpolated  D-function  would  lead  to  singular  integrals  for  the  coefficients 
k  =  I,  k  =  i-I  and  k  =  i~  I.  Therefore  D  will  be  approximated  by  a  polynomial  of  degree 
2  in  the  interval  *•  i  <  i  <  X:+ t.  Splitting  again  into  overlapping  distributions, 
which  this  time  are  parabolic,  and  applying  Cauchy's  pnncipal  value  technique  permits 
integration  of  singular  integrands.  The  coefficient  at  the  middle  of  the  inducing  source 
distribution  is  2iven  bv 


cu  - 


1 


-f.-i 


dD     di 


.-!  .<&  x>  -i 


(4.42) 


with 

<rp_ 
di   '' 


£ \i^±z±._  £_z±zi]_  JLLllL\ -JL—      _Ri 1 

'i-rl-Xi-l    [iTi-Xi-i  Xi+i-Xil         Xi+l-Xl-l   [x  i+i-Xi     ""jf.-A-lJ 
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_  Llfx.^-x.       x. -,._,] 
'- 1 1  —  .  ~— ^- ^—  —  — I  — 
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The  coefficient  at  the  right  of  the  inducing  source  distribution  is  given  by 


"D-iJJ,.1   ft  Jr.  —  ^ 


(4.43) 


with 


2>,--i 


with     -—  = 


di 


D,-i        J-,^1-7,  2(£-x,)        D,_i 


/,-/,_!.    /,  +  !-/,_!  /j  +  l  — /j_i    A-A_i 


for      /,_2  <  $  <  j-.-i 
for     x,-i  <  £  <  j,+1 


,-.-i[ 


m 


J",  —  /,_•> 


J.'-l-l  "J", 


In 


■r,-J\-i 


The  coefficient  at  the  left  of  the  inducing  source  distribution  is  given  by 


•-f.-j 


•  i.i-t-i  — 


dD      c/£ 


*A+i ./,.._,     di  x,  -  £ 


(4.44) 


with 


•  ."     dD 
with     — 
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(        D.-+,        /,-/,_!  2(6-/,)        D, 
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"A-  -f.xi 


la 


■f'l-l  -  x>  )i  ■f'+l  —  ''-1  )  I  ■f|'"~-f(+l  I         J\+j— /,  J-,+2— /,+l  \Xi  —  Xi+2 


As  indicated  above,  the  overall  solution  is  approached  in  an  iterative  process,  in 
which  alternately  viscous  and  inviscid  flow  equations  are  being  solved: 

1.  Calculate  the  external  velocity  distribution  uej  in  an  inviscid  flow  field  by  means 
of  the  conformal  mapping  method.  To  account  for  the  airfoil-thickening  due  to 
viscous  displacement,  specify  a  nonzero  normal  velocity  (blowing  velocity)  at 
the  airfoil  surface. 

2.  March  through  the  boundary  layers  of  both  streamwise  upper  and  lower  surface 
using  the  interaction  law  as  outer  boundary  condition. 

3.  Check  convergence  and  quit  if  the  criterion  satisfied. 

4.  If  the  convergence  criterion  is  not  met.  prepare  for  another  cycle.  Update  the 
product-term  "external  velocity  times  displacement  thickness"  on  the  basis  of 
the  deviation  between  inviscid  and  viscous  external  velocitv  distributions 
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(«e<5-)'-A+l  *(«,«•)<•*     l  +  w( 


H.A 


"e/l 


i. A 


-1 


(4.45) 


where  A.  denotes  the-counter  of  global  iterations.  Further,  compute  the  blowing  velocity 
(  *'vv7r=  waD  transpiration  velocity  )  distribution,  which  serves  as  boundary  condition 
for  the  inviscid  flow  solution 


la.,r)     =   -_(|ljri«.A+l 


(4.46) 


and  proceed  with  the  first  step. 

D.      TURBULENCE  MODELLING 

The  precence  of  additional  unknown  shear  stresses  in  turbulent  flows  requires 
modelling  assumptions  to  balance  the  number  of  unknowns  and  equations.  Eddy 
viscosity  models,  one  of  which  is  used  in  the  present  method.  [Ref.  1?.]  relate  turbulent 
shear  stresses  to  mean  How  quantities  on  an  empirical  basis.  They  draw  their  versatility 
from  the_convemence  of  maintaining  the  same  approach  and  numerical  formulation  for 
both  laminar  and  turbulent  Hows.  According  to  this  formulation  for  wall  boundary 
layer  flows,  the  turbulent  kinematic  eddy  viscosity  is  defined  by  two  separate  formulas. 
one  for  the  inner  region  being  based  on  Van  Dnest's  approach,  and  the  other  for  the 
outer  region  being  based  on  a  velocity  defect  approach 


v*  =   < 


r  ,      7. n  2  ,dn. 

"#[i-«p(-3)]J  1^* 


a  I      [ue- 


'0  dij  it,  t 


for     0  <  ij  <  Uc 
for     ijc  <  y  <  6 


(4.47) 


where 


A  =  2Quf[u — ) 

/    V    ti'Jl  max 


and     7  = 


i  -  o.o\j/d)r' 


Continuity  of  the  turbulent  kinematic  viscosity  is  established  by  defining  >•  as  the 
distance  from  the  wall  where  expressions  for  inner  and  outer  region  do  agree.  Since 
transition  is  not  an  instantaneous  process,  an  intermediate  status  of  flow  is  assumed 
between  the  laminar  and  fully  turbulent  regions.  This  region  is  taken  into  account  by 
introducing  an  intermittency  factor,  which  smears  out  the  step-shaped  change  from 
kinematic  to  eddv  viscosity.   The  formula  here  is  succested  bv  Chen  and  Thvson. 


o 


7tr  =  1  -  exp 


(4.48] 


where   £ 


xtr  denotes   the   transitional   Reynolds   Number  (    ujx;\)tr    i.e., 


Revnolds 


Number  based  on  external  velocity  and  streamwise  location  x[r  at  onset  of  transition. 
A  value  of  1200  was  originally  assigned  to  the  empirical  constant  G,..  .  However, 
numerical  experiments  seem  to  indicate  that  low  Reynolds  number  flows  can  be  better 
modelled  by  values  below  1200  (further  discussion  in  the  next  section).  The  parameter 
a  in  the  outer  resion  formula  is  obtained  from 


0.0168 
a  =  -jzr  =  0.0168y 


1-3 


du/dx]i-* 


du/dy 


where  the  nondimensional  factor  F  denotes  the  ratio  of  total  turbulence  energy 
production  to  the  shear-stress-related  turbulence  energy  production,  evaluated  at  the 
location  of  maximum  shear  stress 


F=  1- 


--,/* 


•du/d. 


a')?'  dnjdij 


(4.49) 


(— a'if'Jmajc 


The  ratio  of  the  time-averaged  quantities  above  can  be  approximated  by  a  function  of 

Rj  —  "<ul\  ~  u'l'lmax 


,?  = 


i2  /2 


J(-u;t»'] 


l  +  2RT{2-RT) 
1  +  Rt 


Rt 


for    RT  <  1 
for    RT  >  1 


(4.50) 


Since  this  turbulence  model  is  not  validated  for  separated  flow,  eddy  viscosities  in 
regions  of  backilow  correspond  to  those  of  adjacent  upstream  station  at  the  same  in- 
coordinate. Transition,  if  not  available  from  other  sources,  currently  is  predicted  by  an 
empirical  data  correlation  expressed  in  terms  of  the  Reynolds  Numbers  based  on 
momentum  thickness  and  streamwise  coordinate  at  onset  of  transition 
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E.       DISCUSSION  OF  THE  COMPLTER  PROGRAM 
1.  Inviscid  Flon  Method 

The  Inviscid  (low  method  used  in  the  interaction  method  is  based  on  conformal 
mapping  which  can  be  divided  into  the  transformation  of  the  region  outside  the  airfoil 
to  the  region  outside  a  unit  circle  and  to  the  solution  of  the  equations  in  the 
transformed  plane.  The  transformation  was  achieved  in  three  parts  which  together 
represent  the  major  computational  effort. 

In  the  first  mapping,  the  airfoil  is  perturbed  slightly  to  make  the  upper  and 
lower  surface  trailing-edge  points  coincide.  This  is  accomplished  using  a  logarithmic 
mapping  function  and  is  necessary  only  in  those  cases  in  which  the  airfoil  trailing  edge 
has  non-zerc  thickness. 

In  the  second  mapping,  the  trailing-edge  corner  is  analytically  removed  by 
applying  the  Karman-Trefftz  mapping. 

In  the  last  mapnng.  the  resulting  quasi-circular  .shape  is  mapped  to  a  perfect 
circle  using  an  iterated  sequence  of  applications  of  the  fast  Founer  transform 
algorithm.  The  calculation  of  the  flow  in  the  transformed  plane  also  makes  use  of 
Fourier  analysis  technique. 

The  major  computational  effort  required  in  the  inviscid  flow  method  is  due  to 
the,  transformation.  It  is  necessary  to  compute  the  transformation  only  once  in  the 
viscous  inviscid  flow  interactions,  so  that  the  computational  expense  due  to  inviscid 
calculations  can  be  held  to  a  minimum.  When  the  angle  of  attack  increases,  care  must 
be  taken  since  the  displacement  thickness  can  become  fairly  large,  approaching  10%  of 
the  airfoil  chord.  In  this  case,  use  of  the  blowing  velocity  (  equation  4.33  )  on  the 
airfoil  surface  produces  a  dividing  streamline  from  the  leading-edge  stagnation  point 
which  approximates  the  edge  of  the  boundary-layer  displacement  thickness.  The 
inviscid  flow  outside  this  dividing  streamline  is  therefore  the  same  as  that  past  the  solid 
body  defined  by  this  streamline.  However,  inside  this  dividing  streamline,  the  inviscid 
flow  is  fictitious.  In  particular,  the  assumption  that  there  is  no  pressure  variation 
across  this  ficitious  region  becomes  invalid  as  the  magnitude  of  blowing  velocity  K 
increases.  The  approach  adopted  here  is.  therefore,  to  evaluate  the  velocity  distribution 
directly  on  the  displacement  surface  while  still  applying  the  blowing  velocity  on  the 
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original  airfoil  surface.  To  avoid  a  discontinuity  of  the  velocity  between  upper  and 
lower  surface,  the  Kutta  condition  is  applied  on  the  displacement  surface.  Requiring 
equal  off-body  pressures  for  the  upper  and  lower  trailing  edge  points,  a  quadratic 
equation  can  be  solved  for  the  unknown  circulation. 
2.   Interactive  Viscous  Flow  Method 

In  this  method,  the  boundary  layer  is  transformed  into  Falkner-Skan  and 
subsequently  semi-transformed  coordinates.  The  numerical  solution  of  the  transformed 
equations  is  performed  using  an  implicit  finite-difference  technique,  which  is  first-order 
accurate  in  the  streamwise  direction  and  second  order  accurate  in  the  normal  direction. 
In  principle,  the  computer  program  calculates: 

1.  The  inviscid  pressure  distribution  for  a  specified  angle  of  attack 

2.  The  boundary  layer,  inclading  the  velocity  profile,  skin  friction,  momentum  and 
displacement  thickness  either: 

1.  subject  to  a  prescribed  pressure  distribution,  or 

2.  subject  to  an  interactive  edge  boundary  condition.  In  this  case,  the  external 
velocity  will  be  a  part  of'the  boundary  layer  results. 

These  calculations  are  performed  for  a  specified  number  of  x-stations  on  the  airfoil  and 
in  the  wake,  a  number  of  sweeps  is  made  on  the  airfoil  in  order  to  obtain  a  converged 
solution.  The  Cebeci-Smith  two  layer  model  is  used  here  for  computing  the  eddy 
viscosity  vf  and  the  transitional  How  region  is  modelled  using-  equation  4.4S.  The 
method  also  employs  a  semi-empirical  formula  to  predict  the  onset  of  transition.  This 
criterion,  equation  4.51,  was  proposed  by  Michel.  The  code  offers  two  possibilities  how 
to  deal  with  transition: 

1.  The  loci  of  transition  are  calculated  by  the  code.  In  this  case  Michel's  criterion 
is  employed  to  predict  the  onset  of  transition.  When  laminar  separation  occurs 
upstream  of  the  calculated  point  of  transition,  then  Michel's  Criterion  is 
disregarded  and  the  onset  of  transition  is  redefined  at  the  point  of  laminar 
separation. 

2.  The  begin  of  transition  is  specified  by  the  user  (fixed).  The  code  has  been 
modified  to  allow  the  onset  of  transition  to  be  within  or  downstream  of  the 
laminar  separation  bubble.  The  previous  version  of  the  code  always  redefined 
transition  at  the  point  of  laminar  separation. 
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F:.2ure  4.5     Lift  curves  of  the  Wortrnann  FX  60-126  airfoil  at  Re=  700.000. 
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(source  of  experimental  results:  Ref.  IS). 


F.       DISCUSSION  OF  THE  RESULTS 

1.  High  Reynolds  Number  Flows 

The  computer  program  was  first  applied  to  high  Reynolds-Number  flows  about  the 
Wortmann  airfoil  FX  60-126.  Figure  4.5  shows  results  for  Reynolds  Numbers  from 
700.000  up  to  2000.000.  It  is  seen  that  the  lift  predictions  are  in  quite  good  agreement 
with  experimental  data. 

2.  Low  Reynolds  Number  Flows 

The  subject  of  low  Reynolds  Number  flows  is  important  to  a  number  of 
practical  problems,  including  remotely  piloted  vehicles,  wind  turbines  and  propellers. 
sail  planes,  human  powered  vehicles,  etc.  Many  boundary  layer  phenomena  which  have 
eluded  accurate  analytical  prediction  occur  within  this  flow  regime  and  are  associated 
with  laminar  separation  and  subsequent  transition  of  a  laminar  free  shear  layer.  Figure 
4.6  shows  the  phenomenological  features  of  the  boundary  layer  on  a  low  Reynolds 
number  airfoil. 

When  the  Reynolds  Number  based  on  momentum  thickness.  Rq.  is  sufficiently  low.  the 
boundary  layer  remains  laminar  up  to,  including,  the  point  of  minimum  pressure  or 
maximum  suction.  At  some  location  between  the  minimum  pressure  and  the 
theoretical  point  of  laminar  separation  the  Reynolds  Number  of  the  boundary  layer 
attains  a  critical  value.  This  is  referred  to  as  the  transition  Reynolds  Number  based  on 
momentum  thickness.  Rotr 

The  provision  of  reliable  information  of  the  flow  characteristics  of  low 
Reynolds  Number  airfoils  'is  hampered  by  the  sensitivity  of  the  flows,  and  particularly 
those  involving  separation  and  transition. 

In  recent  years,  several  theoretical  studies  have  led  to  the  development  o[" 
methods  for  predicting  airfoil  characteristics  at  low  Reynolds  Numbers,  but  there  are 
no  universal  methods  which  can  accurately  predict  and  account  for  a  separation  bubble 
in  the  design  of  efficient  low  Reynolds  Number  airfoils. 

The  viscous,  inviscid  interaction  method  was  applied  to  the  NACA  65-213 

airfoil  at  a  Reynolds  number  of  240,000.    It  was  found  that  the  calculation  would  fail 

to  converge  if  transition  was  predicted  by  Michel's  criterion  (equation  4.51)  and  if  the 

empirical  constant  Gytr  was  chosen  to  be  1200.    Therefore  it  was  decided  to  investigate 

the  influence  of  the  start  of  the  transition  and  of  the  constant  G./tr  on  the  results  bv 

fir 

systematically  varying  both  parameters.  Table  1  of  Appendix  C  shows  the  predicted 
lengths  of  the  separation  bubble  which  is  obtained  if  transition  is  chosen  to  start  at  64, 
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Fieure  4.6     Phenomenoloeical  features  of  the  boundary  laver 
on  the  low  Reynolds  Number  airioil. 


81 


7.4 


u 

1.0 

08 

06 

0£ 

02 

Predicted  Velocity  (IS) 
Measured   Velocity 
i 

i 
i 


I 


•  laminar 


'Seccrcted—~turbr: 
Sucole 


025 


050 


0) 


x/l  '     too 


£sts  b> 


Ficure  4.7     Experimental  result  for  the  NACA  65-213  (bv  Hohcisel  et  al.) 
Ma*-  o.l  .  Rc-=  240.000  .  Tuj-  U.2%  .  AOA  -0.0  deg 

(a).  Velocitv  distribution  and  (b).  I  low  visualization. 
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65.  66.  67,  68,  69.  70,  72,  74.  or  76  %  of  chord  and  if  Gytr  is  chosen  as  10.  20.  40.  SO. 
or  120.  It  is  seen  that  there  are  significant  changes  in  the  length  of  the  separation 
bubble  depending  on  the  chosen  parameter  combination.  This  effect  is  displayed  more 
clearly  in  Figure  4.11  and  4.12. 

An  increase  in  Gvt  increases  the  transition  length  as  well  as  the  length  of  the 
separation  bubble. 

Hoheisel  et  al.  [Ref.  14]  performed  detailed  laser-  Doppler  velocimetry 
measurements  of  this  airfoil  at  a  Reynolds  number  of  240.000.  Their  results  are  shown 
in  Figures  4.".  4.9.  and  4.10  and  it  was  attempted  to  choose  that  parameter 
combination  which  would  give  the  best  agreement  with  the  experimental  results.  If  the 
begin  of  transition  is  chosen  at  "4  %  of  chord  and  if  Gvtr=  20.  then  the  results  shown 
in  Figures  4.13.  4.14  and  4.15  are  obtained.  It  is  seen  that  the  boundary  layer  profiles, 
displacement  and  momentum  thickness  distributions  upstream  of  the  separation  bubble 
are  in  excellent  agreement,  whereas  considerable  deviations  are  found  in  the  bubble. 
Finally.  Figure  4.16  through  4.20  show  the  computed  boundary  layer  velocity  profiles 
for  different  values  of  Gytr  ranging  from  10  to  50 
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V.  CONCLUSION  AND  RECOMMENDATION 

Cebeci's  viscous  inviscid  interaction  program  was  applied  to  the  analysis  of 
steady  two  dimensional  incompressible  flow  past  a  NACA  65-213  airfoil  at  zero  angle 
of  attack  at  a  Reynolds  number  of  240,000.  Predicted  boundary  layer  characteristics 
were  found  to  be  quite  sensitive  to  the  choice  of  boundary  layer  transition  begin  and 
length.  Good  agreement  with  the  experimental  results  of  Hoheisel  et  al  could  be 
obtained  by  proper  choice  of  both  transition  begin  and  length.  Further  detailed 
measurements  and  calculations  for  other  airfoils  at  low  Reynolds  number  are 
recommended  in  order  to  further  validate  the  predictive  capability  of  the 
viscous  inviscid  interaction  method. 
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APPENDIX  A 
FORTRAN  PROGRAM 

CSNOEXT 

*x**x*7r********************************************************** 

*  * 

*  THIS  PROGRAM  CALCULATE  THE  PRESSURE  DISTRIBUTIONS  IN  THE  * 

*  AIRFOIL  AT  ANY  ANGLE  OF  ATTACK,  BY  USING  PANEL  METHOD  IN  * 

*  2-D,INVISCID,  STEADY  FLOW.  * 

*  IMPLEMENTATION  OF  VORTEX  AND  SOUCRE  DISTRIBUTIONS  ARE  USED  * 

*  ALL  VELOCITIES  AND  LENGTH  ARE  NORMALIZED  BY  FREE  STREAM  * 

*  VELOCITY  AND  CHORD  LENGTH  RESPECTIVELY.  * 

*  * 

*  WRITTEN  BY  :  CAPT . INDAF  PHUTUT  SUBROTO  * 

*  NAVAL  POSTGRADUATE  SCHOOL  * 

*  MONTEREY, CA. ,  OCTOBER  1986  * 

*  * 

******  **********************************  ************************* 
C  '   . 

COMMON  /NODE/  N 

COMMON  /AAA/  AA(100,100) , AAX( 100 , 100 ) , AAY( 100 , 100 ) , AAT( 100 , 100 

COMMON  /BBB/  BB(100,100) ,BBX( 100 , 100 ) , BBY( 100 , 100 ) ,BBT(100,100 

COMMON  /STAR/  BSTAR( 100 ) , CSTAR( 100 ) 

COMMON  /ALBE/  ALPA(IOO),  BETA(IOO) 

COMMON  /VEL/  VI ( 100 ) , VXI ( 100 ) , VYI ( 100)      ° 

COMMON  /XY/  X(100)  .Y(100) ,XB(100) ,YB(100) 

COMMON  /SSS/  S(100) ,SIGMA(100) ,SUMB(100)  ' 

COMMON  /ANGLE/  AN , ANG , TPI ,TH( 100) 

COMMON  /PRESS/  CP(IOO) 
C 

C INPUT  DATA  :  #NODE ,  #  AOA   

C 

READ (5,1)  NN,AN 
1    FORMAT(I10,F10.4) 

PI   =  4.*ATAN(1.0) 

ANG   =  AN/180.*PI 

VXA   =  -COS (ANG) 

VYA   =  -SIN (ANG) 

TPI  =  .5/PI 

DO  10  I   =  1,NN 

READ(5,5)  XB(I),YB(I) 
10     CONTINUE 
5   FORMAT (2F10. 5) 
C 
*********************************************************** 

*  * 

*  '    COMPUTE  THE  ANGLES  AND  LENGHTS  OF  EACH  PANELS,       * 

*  CALCULATE  THE  COORDINATE  OF  CONTROL  POINTS  * 

*  * 

*********************************************************** 
C 

N  =  NN-1 
SS   =  0.0 

DO  15  I  =  1,N 
X(I)    =  .5*(XB(I)  +  XB(I+1 
Y(I)    =  .5*(YB(I)  +  YB(I+1 

TH(I)   =  ATAN2  (YB ( 1+1 ) -YB ( I )  ,  XB( 1+1 ) -XB( I ) ) 
S(I)    =  SQRT((XB(I+1)-XB(I))**2  +  (YB(I+1 )-YB(I) )**2) 
SS      =  SS  +  S(I) 
15  CONTINUE 
C 
*********************************************************** 

*  * 

*  COMPUTE  TIME  INDEPENDENT  INFLUENCE  COEFFICIENTS      * 
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*     FOR  NORMAL,  TANGENTIAL  ,X  AND  Y  DIRECTION  * 


DO  25  1=  1, 

N 

DO  25  J=  1,N 

IF (I  .EQ.  J)  THEN 

AA(T,J)  =   0.5 

BB(I,J)  =   0.0 

AAX  < 

I*J 

i=  -0.5*SIN(TH 

BBX! 

I, J 

i=   0.5*COS(TH 

AAYI 

I, J 

i=   0.5*COS(TH 

BBY( 

I, J 

i=   0.5*SIN(TH 

AATI 

I.J 

i=   0.0 

BBT( 

I.J 

i=   0.5 

ELSE 

A   =  -(X(I)-XB(J))*COS(TH(J) 

B   =    XI 

-XB( 

J))**2  +  (Y(I 

(Y(I)-YB(J))*SIN(TH(J)) 
(J))f^ 


AA  I,J 

BB(] 

:,J 

AAXI 

i. 

EBX< 

i* 

AAY< 

I* 

BBY( 

i, 

AAT( 

i. 

B3T( 
END  IF 

;i, 

25 

CONTINUE 

-YB(J))**2 
C   =   SIN(TH(I)-TH(J 
D   =   COS(TH(I)-TH(J 
E   =   (X(I)-XB(J))*SIN(TH(J))  -  (Y(I)-YB(J))*COS(TH(J)) 
F   =   ALOG(1.0  +  S(J)*(S(J)+2.*A)/B) 
G   =   ATAN2  (E*S(J)  ,  B+A*S(J)) 

TPI*(.5*C*F  -  D*G 
TPI*( .5*D*F  +  C*G 
J)  =  TPI*(-.5*COS(TH(J))*F  +  SIN(TH(J))*G 
J  =  TPI*(-.5*SIN(TH(J))*F  -  COS(TH(J))*G 
J)  =  BBX(I,J) 
J)  =  -AAX(I,J) 
J  =  -BB(I,J' 
J)  =   AA(I,J 


C 

*  SETUP  MATRICES, SET  GAMA=1.0,  SOLVE  THE  SYSTEM       * 

*  BY  USING  GAUSSIAN  ELIMINATION  WITH  PARTIAL  PIVOTING  * 

*  FOR  TWO  RIGHT  HAND  SIDES.  * 

*  * 

c  " 

DO  68  I  =1,N 

SUMB(I)  =0.0 
DO  680  J  =1,N 
680      SUMB(I)    =  SUMB(I)  +  BB(I,J) 
BSTAR(I)   =  -SUMB(I) 

CSTAR(I)   =   -VXA  *  SIN(TH(I))   +  VYA  *  COS(TH(I)) 
AA(I,N+1)  =  BSTAR(I' 
AA(I,N+2)  =  CSTAR(I 
63   CONTINUE 
C 

CALL  GAUSS (2) 
C 

DO  500  I  =1,N 

ALPA(I)  =  AA(I,N+1) 
BETA(I)  =  AA(I,N+2) 
500   CONTINUE 
C 

C   MEET  KUTTA  CONDITIONS  , SOLVE  FOR  VORTEX  STREGNTH(GAMA) 

C       USING  QUADRATIC  EQUATION  

C 

AX1  =  0.0 
AXN  =0.0 
AY1  =0.0 
AYN  =0.0 
BX1  =  0.0 
BXN  =0.0 

99 


BY1  =  0.0 
BYN  =0.0 

DO  510  J  =1,N 

AX1  =  &X1  +  AAX(1,J)*ALPA(J)+BBX(1 
AXN  =  AXN  +  AAX(N,J)*ALPA(J)+BBX  N,J' 
AY1  =  AY1  +  AAY(l,J)*ALPA(J)+B3Y(l,J 
AYN  =  AYN  +  AAY(N,J)*ALPA(J)+BBY(N,j' 
BX1  =  BX1  +  AAX(1,J)*BETA(J' 
BXN  =  BXN  +  AAX(N,J)*BETA( J' 
BY1  =  BY1  +  AAY(1,J)*BETA(J' 
BYN  =  BYN  +  AAY(N,J)*BETA(J, 
510   CONTINUE 
C 

EE  =  AX1**2+AY1**2  -  AXN**2-AYN**2 

PP  =   AX1*BX1+AY1*BY1  -  AXN*BXN-AYN*BYN  +(AXN-AX1 )*VXA 
+        +  (AYN-AY1)*VYA 

QQ  =  BX1**2+BY1**2-BXN**2-BYN**2  +2 . *(BXN-BX1)*VXA  +2.* 
+  (BYN-BY1)*VYA 

C 

R  =   PP*PP  -  EE*QQ 
GAMA1    =  (-PP  +  S0RT(R))/(EE) 
GAMA2   =  (-PP  -  SQRT(R))/(EE) 
IF(ABS(GAMA1)  .GT.  ABS (GANA2) )  THEN 
GAI-IA    =  GAMA2 
ELSE 

GAMA    =  GAMA1 
END  IF 
C 

C     SOLVE  THE  SOURCE  STRENGTH   SIGMA(J)  

C 

DO  550  J  =1,N 

SIGMA(J)=  GAMA*ALPA(J)  +  3ETA(J) 
■  550   -  CONTINUE 
C 

CALL  CPRESS(GAMA) 
C 
C 

C   PRINT  LIFT  COEFFICIENT  AND  MOMENT  COEFFICIENT  ABOUT 

C  LEADING  EDGE  

C 

CFX  =0.0 
<   CFY   =0.0 
CM    =0.0  • 
DO  110  I  =  1,N 

DX  =  XB(I+1)  -  XB(I' 
DY  =  YB(I+1)  -  YB(I' 
CFX=  CFX  +  CP(I)*DY 
CFY=  CFY  -  CP(I  *DX 
CM  =  CM  +  CP(I)*(DX*X(I)+DY*Y(I)) 
110  CONTINUE 

CL    =  CFY*COS(AMG)  -  CFX*SIM(ANG) 
CD    =  CFX*COS(ANG)  -  CFY*SIN(ANG) 
WRITE (8, 120)  CL 
WRITE (8 ,125)  CD 
WRITE (3 ,130)  CM 
120  FORMAT (///, 1 5X, ' CL      =  \F10.5) 
125  FORMAT(/,15X,   'CD      =',F10.5) 
130  FORHAT(/,15X,   ' CMLE     =',F10.S) 
PRINT*, ' "COMPUTATION  COMPLETED" ' 
STOP 
END 

************************************************************ 

*  * 

*  SUBROUTINE  GAUSS  * 

*  SOLVE  SIMULTANEOUS  EQUATION  WITH  TWO  RIGHT  HAND  SIDES  * 

*  BY  GAUSS  ELIMINATION  WITH  PARTIAL  PIVOTING  * 

*  SOLUTIONS  STORES  IN  COLUMNS  NEQMS+1  AND  NEQNS+2  OF  * 

*  MATRIX  AA(NEQNS, NEQNS  +  2)  .  * 

*  (AA)   =  COEFFICIENT  OF  AUGMENTED  MATRIX  * 
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*  NEONS  =  NUMBER  OF  EQUATIONS  * 

*  KRHS   =  NUMBER  OF  RIGHT  HAND  SIDES  * 

*  * 

SUBROUTINE  GAUSS (NRHS) 
C 

COMMON  /NODE/  NEQNS 

COMMON  /AAA/  AA(100 , 100) , AAX(100 , 100) , AAY(100 , 100) ,AAT(100 , 100) 

N?         -   NEONS +1 

::TOT       =  NEQNS+NRKS 

c 

C  GAUSS  REDUCTION 

DO  40   I   =  2, NEQNS 
C 

C  SEARCH  FOR  LARGEST  ENTRY  IN  (I-l)TH  COLUMN 

C  ON  OR  3EL0W  MAIN  DIAGONAL 

C 

IM         =  1-1 

IMAX       =  IM 

AMAX       =  ABS(AA(IM,IM)) 

DO  10    J  =  I, NEQNS 

IF (AMAX  .GE.  ABS(AA(J,IM)))   GO  TO  10 
IMAX    =  J 

AMAX    =  ABS(AA(J,IM)) 
10   CONTINUE 
C 

C  SWITCH  (I-l)TH  AND  IMAXTH  EQUATIONS 

C 

c 

IF(IMAX  .NE.  IM)   GO  TO  30 
DO  20    J  =  IM,NTOT 

TEMP        =  AA(IM,J) 
AA(IM,J)      =  AA(IMAX,J) 
AACMAX,J)    =  TEMP 
20   CONTINUE 
C 

C  ELIMINATE  (I-l)TH  UNKNOWN  FROM  ITH  THRU 

C  (NEQNS )TH  EQUATIONS 

C 

30      DO   40      J   =    I, NEQNS 

R  =   AA(J,IM)/AA(IM,IM) 

DO   40  X  =  T,NTOT 

AA(J,XV  =  AA(J,X)-R*AA(IM,K) 

40      CONTINUE 
C 

C  3ACX   SUBSTITUTION 

C 

DO  70   X  =  NP,NTOT 

AA(  NEONS,  X)   =  AiU  NEONS, K)/AA (NEQNS, NEQNS) 
DO  50 ~    L  =  2, NEONS 

I        =  NEQtiS+1-L 
I?       =  1+1 
DO  50   J  =1?, NEQNS 
50   AA(I,X)     =  AA(I,X)-AAn,J)*AA(J,K) 
60   AA(I,X)     =  AA(I,X)/AA(I,I) 
70   CONTINUE 
RETURN 

z::: 

c 
c 

*  * 

*  SUBROUTINE  CPRESS  * 

*  CALCULATE  PRESSURE  COEFFICIENTS  AND  TOTAL  VELOCITY     * 

*  AT  MID  POINTS  (CONTROL  POINTS)  * 

*  -k 

SUBROUTINE  CPRESS (GAMA) 
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COMMON  /AAA/  AA( 100 , 100) , AAX( 100 , 100 ) , AAY( 100 , 100 ) , AAT( 100 , 100 
COMMON  /BBB/  BB( 100 , 100) ,BBX(100 , 100) ,BBY( 100 , 100) ,BBT( 100 , 100 
COMMON  /PRESS/  CP(IOO) 
COMMON  /NODE/  N 

COMMON  /ANGLE/  AN, ANG,TPI ,TH( 100) 
COMMON  /SSS/  S(100) . SIGMA( 100 ) , SUMB( 100 ) 
COMMON  /VEL/  VI ( 100 ) , VXI ( 100 ) , VYI ( 100 ) 
COMMON  /XY/  X(1Q0),Y(100),XB(100),YB(100) 
WRITE (S, 5) 

5  FORMAT(//,15X,  '  NACA  23012'  ) 

WRITE (8, 5)  N 
WRITE (3, 7)  AN 

6  FORMAT (15X, 'NUMBER  OF  PANELS  =',I3) 

7  FORMAT (15X,  'ANGLE  OF  ATTACK   =',F7.4) 

WRITE (8, 20) 
DO  12  I  =  1,N   ■ 
SUMM   =0.0 
DO  13  J  =  1,N 

SUMM  =  SUMM  +  AAT(I,J)*SIGMA(J)+GAMA*BBT(I,J) 
13   CONTINUE 

VI (I)  =  SUMM  +  COS(TH(I)-ANG) 
CP(I)   =  1.0  -  VI(I)**2 
•  WRITE(3,30)  I,X(I)/Y(I)/VI(I),GAMA,SIGMA(I)/CP(I) 
WRITE(6,31)  X(I),-CP(I) 
12   CONTINUE 
20   FORMAT(///3X, 'PNL(I) ' , 3X, ' X(I ) ' , 7X, ' Y(I ) ' , 6X, ' VEL (I ) ■ , 5X, 


1  'GAMMA'  ,5X, 'SIGMA(I) '  ,4X,  ' CP ( I )  '  /J) 

30  FORMAT(2X,I3,5X/F3.5/3X/F8.5,3X,F8.5,3X,F8.5,3X/F8.5,3X,F8.5) 

31  FORMAT (2F1 0.5) 

RETURN 
.END 


102 


APPENDIX  B 
PROGRAM  OUTPUT 


PML(I)    X(I) 


SOURCE  AND  VORTEX  PANEL  SOLUTION 

NACA  23012 

NUMBER  OF  PANELS  =  34 

ANGLE  OF  ATTACK   =12.0000 


Y(I) 


VEL(I) 


GAMMA 


SIGMA(I) 


CP(I) 


1 

0.97500 

-0.00350 

-0.89112 

0.39035 

2.28747 

0.20590 

2 

0.92500 

-0.00965 

-0.88365 

0.39035 

2.31831 

0.21917 

3 

0.85000 

-0.01695 

-0.88620 

0.39035 

2.40999 

0.21464 

4 

0.75000 

-0.02580 

-0.88782 

0.39035 

2.41656 

0.21177 

5 

0.65000 

-0.03335 

-0.87216 

0.39035 

2.42538 

0.23933 

6 

0.55000 

-0.03920 

-0.85087 

0.39035 

2.43692 

0.27602 

7 

0.45000 

-0.04325 

-0.82394 

0.. 39035 

2.44450 

0.32113 

8 

0.35000 

-0.04470 

-0.77104 

0.39035 

2.47951 

0.40549 

9 

0.27500 

-0.04370 

-0.71339 

0.39035 

2.50195 

0.49107 

10 

0.22500 

-0.04125 

-0.66069 

0.39035 

2.54546 

0.56343 

11 

0.17500 

-0.03735 

-0.57694 

0.39035 

2.60630 

0.66714 

12 

0.12500 

-0.03210 

-0.46811 

0.39035 

2.65463 

0.78087 

13 

0.08750 

-0.02765 

-0.36096 

0.3903.5 

2.64815 

0.86971 

14 

0.06250 

-0.02435 

-0.25392 

0.39035 

2.62380 

0.93553 

15 

0.03750 

-0.01985 

-0.00657 

0.39035 

2.65274 

0.99996 

16 

0.01875 

-0.01470 

0.41263 

0.39035 

2.58301 

0.82974 

17 

0.00625 

-0.00615 

1.41005 

0.39035 

.  2.6722*4 

-0.98825 

13 

0.00625 

0.01335 

2.61069 

0.39035 

-0.18981 

-5.31563 

19 

0.01875 

0.03140 

2.59588 

0.39035 

-1.31312 

-5.73860 

20 

0:03750 

0.04260 

2.32749 

0.39035 

-1.62967 

-4.41720 

21 

0.06250 

0.05355 

2.17813 

0.39035 

-1.84536 

-3.74423 

22 

0.08750 

0.06115 

2.07578 

0.39035 

-2.00811 

-3.30837 

23 

0.12500 

0.06810 

1.90685 

0.39035 

-2.23615 

-2.63608 

24 

0.17500 

0.07345 

1.74625 

0.39035 

-2.40177 

-2.04941 

25 

0.22500 

0.07550 

1.63801 

0.39035 

-2.46583 

-1.63307 

26 

0.27500 

0.07575 

1.55856 

0.39035 

"2.50044 

-1.42910 

27 

0.35000 

0.07345 

1.45669 

0.39035 

-2.55339 

-1.12194 

23 

0.45000 

0.06775 

1.35494 

0.39035 

-2.58751 

-0.83586 

29 

0.55000 

0.05940 

1.27629 

0.39035 

-2.61244 

-0.62891 

30 

0.65000 

0.04915 

1.20826 

0.39035 

-2.62994 

-0.45989 

31 

0.75000 

0.03720 

1.13842 

0.39035 

-2.65075 

-0.29600 

32 

0.35000 

0.02380 

1.07219 

0.39035 

-2.65258 

-0.14953 

33 

0.92500 

0.01300 

1.01371 

0.39035 

-2.57577 

-0.02761 

34 

0.97500 

0.00460 

0.89101 

0.39035 

-2.59479 

.  0.20611 

CL 

1 

.57532 

CD 

=   -0 

.63272 

CMLE     =   -0, 

.41146 
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APPENDIX  C 
TABLE  I 

TABLE  2 
EFFECT  OF  GGTR  AND  XTRL  ON  THE  BUBBLE  LENGTH 


GGTR=10 

20 

40 

SO 

120 

XTRL  =  .64 

.6567-.6745 

.6567-.6920 

.6567-.7093 

.6567-J263 

.6387-7593 

XTRL  =65 

.6567-. 6920 

.6745-."093 

.6567-.7263 

.6387-. 7429 

.6387-.7752 

XTRL  =.66 

.6745-7093 

.6567-.7093 

.6567-.7263 

.6387-.7593 

.6204-.7908 

XTRL  =  .67 

.6745-. 7093 

.6745-.7093 

.63S7-.7429 

.6204-. 7752 

.6204-.8207 

XTRL  =  .68 

- 

.6567-.7263 

.6387-. 7593 

.63S7-.7593 

.6204-8060 

.6020-.8489 

XTRL  =  .69 

.63S7-.7429 

.6387-."593 

.6204-.7752 

.6020-.8350 

.5834-. 8876 

XTRL  =.70 

.63S7-.7593 

.6204-.7752 

.6204-.8060 

.  .6020-.8623 

.6204-.8060 

XTRL  =  .72 

.6204-.7908 

.6204-.8060 

.6020-.8350 

.5834-.S752 

.5647-.9216 

.9415-.9740 

(turb) 

XTRL  =.74 

.6020-.8060 

.6020-.8350 

.5834-.S623 

.5647-.9319 

. 

XTRL  =.76 

.6020-.8350 

.5834-.8623 

.5647-.S752 

- 

- 

NACA  65-213 

Revnolds  No. 

AOA 

XTRL 

XTRL 

GGTR 


=  240.000 

=  0.0  deg 

=  Begin  of  transition  (  Lpper  Surface  ) 

=   Begin  of  transition  (  Lower  Surface  fixed  at  0.4  ) 

=  Empirical  Constant  (G-y    ) 
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APPENDIX  D 
TABLE  II 

TABLE  3 

EFFECT  OF  GGTR  AND  XTRL  ON  THE  SHAPE  FACTOR(H) 
AT  POINT  OF  ZERO  SKIN  FRICllON 


GGTR=10 

20 

40 

80 

120 

XTRL  =  .64 

2.830 

2."80 

2.726 

2.7423 

2.722 

XTRL  =.65 

2."262 

2.-45 

2.684 

2."024 

2.-013 

XTRL  =  .66 

• 

2.586 

2.624 

2.660 

2.691 

2.682 

XTRL  =.67 

2.b04 

2.6S7 

2.5941 

2.653 

2.78 

XTRL  =.68 

2.95 

2.478 

2.39 

2.6478 

2.589 

XTRL  =  .69 

2.8" 

2.776 

2.644 

2.5554 

2.5844 

JCTR.U-.7Q 

2  5524 

2.523 

2.534 

2.5152 

2.478 

XTRL  =."2 

2.54 

2.583 

2.60 

2.487 

2.480 

XTRL  =."4 

• 

2.575 

2.609 

2.692 

2.669 

• 

XTRL  =  ."6 

2.624 

2.667 

2.82 

NACA  65-213 

Reynolds  No. 

AOA 

H 

XTRL 

XTRL 

GGTR 


=  240.000 
=  0.Q,de2 

=  6   0" 

=   Begin  of  transition  (Lpper  Surface) 

=  begin  of  transition  (Lower  Surface  fixed  at  0.4) 

=  Empirical  Constant  (G„,    ) 
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APPENDIX  E 
TABLE  III 

TABLE  4 
EFFECT  OF  GGTR  AND  XTRL  ON  THE  DRAG  COEFFICIENT  (CD) 


GGTR=10 

20 

40 

80 

120 

XTRL  = 

.64 

TE:.01025 
WK.-.01064 

.01015 
.01054 

.010029 
.010423 

.009876 
.010274 

.009813 
.010225 

XTRL  = 

.65 

TE:.01022 
\VK:.0I062 

.01013 
.01052 

.010022 
.010416 

.009900 
.010301 

.009785 

010193 

XTRL  = 

.66 

TE:.01010 
\VK:.01049 

.010112 
.01050 

.010018 
.010398 

.009845 
.010248 

.009741 
.010149 

' XTRL  = 

.67 

TE:.01016 
WK:.0l055 

.01010 
.01050 

.00926 
.010316 

.009785 
.010185 

.009816 
.010236 

XTRL  = 

.68 

TE:.01<;16 
\VK:.01055 

.009980 
.010383 

.009944 
.010338 

.009846 
.010246 

.009865 
.010312 

XTRL  = 

.69 

TE:.01001 
WK:.01050 

.010015 
.010405 

.009948 
.010339 

.009886 
.010316 

.010014 

.010467 

XTRL  = 

.70 

TE:.010ll 
WK:.01049 

.010015 
.010403 

.009900 
.010334 

.009964 
.010382 

.010173 
.010729 

XTRL  = 

.72 

TE:.01008 
VVK:.0l047 

.00997 
.010352 

.009981 
.010380 

.010034 
.010481 

.010644 
.011202 

XTRL  = 

.74 

TE:.01004 
WK:01042 

.01005 
.010405 

.010017 
.010411 

.010489 
.010947 

- 

XTRL  = 

.76 

TE:.00999 
WK:.0l037 

.010024 
.010434 

.010222 
.010672 

. 

- 

NACA  65-213 

Revnolds  No. 

AOA 

TE 

VVK 

XTRU 

XTRL 

GGTR 


=  240.000 

=  0.0deg 

=  Cd  at  the  trailing  edge 

=  Cd  at  the  wake 

=   Begin  of  transition  (Lpper  Surface) 

=   Begin  of  transition  (Lower  Surface  fixed  at  0.4) 

=  Empirical  Constant  (G^    ) 
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APPENDIX  F 
TABLE  IV 

TABLE  5 
EFFECT  OF  GGTR  AND  XTRL  ON  THE  LIFT  COEFFICIENT  (CL) 


GGTR=10 

20 

40 

80 

120 

XTRL-.64 

.1550 

.1555 

.1559 

.1562 

.1560 

XTRL  =.65 

.1557 

.15602 

.1565 

.1569 

.1571 

XTRL  =\66 

.157] 

.15630 

.1573 

•       .1581 

.1589 

XTRL  =6" 

.1568- 

.15692 

.1588 

.1598 

.1599 

XTRL  =68 

1  570 

.1588 

.1594 

.1609 

.1630 

XTRL  =  .69 

.1580 

.15963 

.1609 

.1632 

.1693 

CTRL' -.70 

.1596 

.16072 

.1625 

.1677 

.1800 

XTRL-.72 

.1624 

.16464 

.1687 

.1815 

.2074 

XTRL  =."4 

.1671 

.1708 

.1784 

.2045 

- 

XTRL  =."6 

.1740 

.1805 

.1910 

NACA  65-213 

Revnolds  No. 

AOA 

XTRL 

XTRL 

GGTR 


=  240.000 

=  0.0  deg 

=  Begin  of  transition  (Lpper  Surface) 

=   Begin  of  transition  (Lower  Surface  fixed  at  0.4) 

=  Empirical  Constant  (G„,    ) 
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